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ABSTRACT 


Passive microwave satellite data provides all-weather, continuous 
observations of the polar ice cover. Multichannel radiometer data sets extending 
back to 1978 are readily available and distributed on CD-ROM media. Current 
algorithms can extract information on ice cover, concentration and type. Two 
applications of these data sets for polar studies are explored: 1) the possible 
detection of altered ice morphology due to anomalous tidal currents in the 
Yermak Plateau region; and 2) the tracking of multiyear ice features to infer ice 
motion. 

No evidence of altered ice morphology in the Yermak Plateau region was 
detected. Either the reported anomalous currents do not affect ice morphology or 
the limitations inherent to passive microwave imagery (poor spatial resolution) 
prevent detection. On the other hand, distinct multiyear ice features are observed 
to move in the direction of mean ice flow and are well correlated with the motion 
of nearby data buoys. Short term ice motion vectors are unreliable because of the 
low spatial resolution of the passive microwave data. However, long term 
averages are indicative of ice flow trajectories and can be used, for example, to 
increase our knowledge of the interannual variability of heat and salinity budgets 
in specific regions, including Fram Strait, the principal outflow region of the arctic 


basin. 
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I. INTRODUCTION 


All-weather, continuous observations of the arctic ice cover are important for 
commercial and naval operations and for scientific studies such as global climate 
change. However, such observations are generally difficult to obtain. In-situ 
measurements from ice camps and ships are expensive to acquire and spatially 
limited. Aircraft measurements can expand the areal coverage, but asset and 
resource limitations prevent continuous large-scale monitoring of the ice pack. 
Satellite platforms are the only practical means of providing continuous coverage. 
Of the passive sensing methods, only microwave emissions penetrate the clouds 
and darkness that shield much of the arctic ice from visual and infrared 
observations. Passive satellite microwave data are operationally available from 
Special Sensor Microwave/Imager (SSM/I) sensors on the Defense 
Meteorological Satellite Program (DMSP) series satellites (Hollinger et al., 1987). 
Historical data are also available to researchers on CD-ROM distributed by the 
National Snow and Ice Data Center (NSIDC). 

The objectives of this study are to explore the use of passive microwave 
Satellite measurements for arctic sea ice studies. Two applications are examined: 

¢« The detection of altered ice morphology due to anomalous 
diurnal tidal currents observed over the Yermak Plateau. If 
successful, this technique can be applied in other polar regions to 
expand our knowledge of the local tidal structure and its effect on 
sea ice cover. 

¢ Determination of ice flow by tracking the movement of floes 


exhibiting different morphologies from the surrounding ice field. 


This can augment the Arctic Buoy Program and increase our 
understanding of ice dynamics. 

Chapter II presents background material on the SSM/I, data sources and sea 
ice microwave radiation characteristics. Determination of sea ice concentrations 
from passive microwave data are also discussed. Chapter III examines the 
anomalous diurnal currents observed over the Yermak Plateau. Previous 
observations and models prove the existence and explain the physics behind 
these currents. The effect of these currents on the ice cover is examined through 
passive microwave emissions. Chapter IV shows how ice features can be tracked 
to infer general ice drift patterns in the arctic. When combined with other data, 
this technique can help researchers understand ice and salinity budgets in the 
polar region. Chapter V summarizes the results and closes with a discussion of 


further research that could be attempted using these data. 


Il. BACKGROUND 


A. THE SSM/I 

The SSM/I is flown on the Defense Meteorological Satellite Program (DMSP) 
block 5D-2 F8 satellite, launched on June 19, 1987. It is in circular sun- 
synchronous near-polar orbit at an altitude of approximately 833 km, an 
inclination of 98.8° and a period of 102 minutes. This produces 14.1 full orbit 
revolutions per day with a 0613 local ascending node equatorial crossing. The 
SSM/I's conical scan covers a 1400 km swath which provides 24 hour near-global 
coverage. The orbit inclination precludes coverage of a 280 km radius region 
centered on each pole extending to about 87.5°N. Equatorial regions require 72 
hours for full coverage, however overlapping orbits provide full daily coverage 
poleward of 45°N. Hollinger et al. (1990) provide a good discussion of SSM/I 
operating characteristics. 

The SSM/I uses four frequencies to measure seven channels of linearly 
polarized radiation as shown in Table 1. Footprint resolution increases with 
increasing frequency. Unfortunately, the high resolution 85 GHz vertically 
polarized (85V) channel was damaged when the satellite overheated in the winter 
of 1987. The channel was useless by the end of January 1989 (Hollinger et al., 
1990). The National Snow and Ice Data Center (NSIDC) recently reported that 
the 85 GHz horizontally polarized (85H) channel is becoming increasingly 
unreliable and is basically useless (Taylor, 1990). Poe and Conway (1990) 


reported the current SSM/I geo-location error to be approximately 10-15 km. 


NSIDC applies a correction algorithm that reduces this error to 8 km (NSIDC, 
WON 


TABLE 1: SSM/I FREQUENCIES 
Ch 1 | Frequency | Wavelength | polarizati Resolution 


Vertical 


Horizontal 33 


ak 
‘ne 





SSM/I data are transmitted to Telemetry and Data Relay Stations on each orbit. 
These data are then relayed via satellite to Navy and Air Force processing centers 
which produce near real time global maps of sea ice location, age, type and other 
environmental parameters. Raw data are also continually broadcast to Air Force 
and Navy tactical units positioned within the orbital swath of the satellite that are 


equipped with SMQ-11 satellite receivers. 


Data record files are created from raw radiometer measurements. Temperature 
data records (TDRs) are antenna temperatures which are calibrated each scan 
using a 300° K black-body radiator and the 3° K cosmic radiation as reference 
temperatures. An antenna pattern correction is applied to the TDRs to correct for 
energy received in the sidelobes, cross-polarization coupling and feedhorn 
spillover loss (NSIDC, 1990). These main beam brightness temperatures are 
placed in sensor data record (SDR) files. Environmental parameters are 
determined by applying appropriate algorithms to the SDRs and are stored in 
environmental data records (EDRs). The primary archiver of SSM/I data in orbital 
format is the NOAA/National Environmental Satellite, Data, and Information 
Service, Satellite Data Services Division (NOAA/NESDIS/SDSD). TDRs, SDRs 


and EDRs are also forwarded to the Naval Research Laboratory for archival. 


B. SEA ICE RADIATION CHARACTERISTICS 

Sea ice is a complex material. A good comprehensive overview of both its 
physical and microwave radiation characteristics can be found in Parkinson et al. 
(1987). The following is a general discussion of these properties. 

1. Physical Properties 

Sea ice covers approximately seven percent of the word's oceans. It 1s 

inhomogeneous, anisotropic and never in equilibrium. Sea ice can simultaneously 
contain all three phases of matter: solid ice, liquid brine pockets and air cavities. 
The fractions of the constituent phases change over time. New ice is thinner and 
saltier than first year ice, which transforms after a melt cycle (summer Season) into 
fresher, thicker and less dense multiyear ice containing numerous air pockets. Ice 


can be covered with a layer of dry or wet snow of varying thickness; summer 


melting forms ponds of liquid water over the ice sheet. Figure 1 shows the 
different physical characteristics of first and multiyear ice. 

The ice composition is initially determined by environmental factors at 
the time it is formed. These factors include air temperature, wind speed, ocean 
currents, near surface salinity and weather (Weeks, 1981). The ice composition 
changes continuously due to weather, ice advection and fracturing or ridging 
effects. Warmer summer temperatures are especially important as melt water 
percolates through the ice and brine pockets are flushed. 

2. Radiative Properties 

The radiative properties of sea ice are determined by the emissivity 
differences between ice and water. Emissivity is the ratio of emittance (J m-2 s-!) 
per unit wavelength to that of a blackbody at the same temperature. The 
Rayleigh-Jeans approximation (linear relationship between emittance and 
temperature) applies at microwave wavelengths. The emissivity (€) of a substance 
can be expressed as: 


Te 
a T sfc (1) 





where Te is the radiated microwave brightness temperature, proportional to the 
emittance, and Tsfc is the physical temperature of the substance. The SSM/I 
radiometer senses TB, the brightness temperature at the satellite which is related 
to Te after attenuation and contamination by atmospheric effects are considered. 
Brightness temperatures can be used to calculate ice concentrations and 
differentiate between open water, first year ice and multiyear ice based on 


emissivity variations of the substances. To a lesser extent, snow and water- 


covered ice can also be observed. Figure 2 summarizes the emissivity differences 
between sea water, first year ice and multiyear ice as a function of polarized 
frequency. 

The ice composition controls the dielectric properties that influence the 
emissivity of the ice. Following Comiso (1983), the principal difference between 
sea ice and freshwater ice are brine inclusions in the former. These brine pockets 
are small, typically 3-5 mm by 0.25 mm (Vant et al., 1978). The brine content of 
sea ice depends on freezing conditions and subsequent temperature changes and 
varies from 16 psu for first year ice to less than 1 psu for multiyear ice. The liquid 
brine effects the dielectric property that determines emissivity and the skin depth, 
which is the depth from which most of the observed radiation emanates. In saline 
ice the skin depth is on the order of the radiation wavelength. For desalinated ice 
this depth can be many wavelengths (Gloersen et al., 1973). 

The fundamental difference between first year and multiyear ice is the 
presence of brine pockets in the former and their replacement in the freeboard 
layer by air in the latter. Because of the small skin depth, microwave emissions 
from the more saline first year ice emanate from the snow/ice interface; emissions 
from fresher multiyear ice emanate from a thicker ice layer and are subject to 
scattering by the air pockets. This scattering reduces the net radiation from the 
ice surface thus lowering the brightness temperature (Comiso, 1983). 

Snow cover also effects microwave radiation. Snow insulates the ice 
from the atmosphere, leading to warmer ice from oceanic heat conduction and 
therefore slightly higher brightness temperatures. This may be compensated by 
the scattering effect of snow ice crystals, especially at higher frequencies. 


Emissions from the crystals are negligible (Comiso et al., 1982). During the 


summer melt season, however, each crystal is surrounded by a thin layer of water. 
This wet snow cover becomes the primary radiation source, significantly lowering 
brightness temperatures due to water's low emissivity. This effect makes it nearly 
impossible to discriminate wet snow-covered multiyear ice from first year ice in 
the summer. 

Ice thickness also influences brightness temperatures. At microwave 
frequencies, thin ice (radiation wavelength or less) is contaminated by the coldly 
radiating water below, lowering the brightness temperature. This may be 
compensated by the higher salinity and warmer physical temperature of the thin 
ice versus the relatively fresher and colder thick ice. 

Melt ponds are formed by solar heating of the ice cover and are 
characterized by thin layers of water scattered among solid ice regions. The size 
and depth of the ponds depend on the incident solar flux, ice salinity and ice 
roughness (Goroch and Fett, 1990) If the thin water layer is thicker than the skin 
depth of the radiation (typically a few mm to one cm), the radiation 1s emitted 
primarily by the liquid and has a much colder brightness temperature than the ice. 
For thinner layers the underlying ice will contribute some radiation, raising the 
brightness temperature. The variation of brightness temperature depends on the 
layer depth and the emission frequency. For depths on the order of a few 
millimeters the 19 and 22 GHz channels will be unaffected (Goroch and Fett, 
1990). As the frequency increases, contamination by the water layer also 
increases. The total effect on brightness temperature is a function of areal extent 
within the radiometer footprint and melt pond depth. 

Other factors affecting microwave emissivities include the physical 


temperature of ice and undulations caused by ridging and refreezing of melt 


ponds. As with liquid melt ponds, these conditions vary spatially and temporally 
and could cause observable changes in ice emissivity. 
3. Atmospheric Effects 

Important atmospheric constituents which affect microwave radiation 
are water vapor, oxygen and liquid water content. Seasonal climatological 
constants developed from physical models or observations are generally used in 
algorithms to quantify their cumulative effect on brightness temperatures. 

Interaction with non-precipitating clouds in the arctic region is 
negligible due to the small size of the droplets compared with the microwave 
wavelengths (Comiso, 1983). Liquid precipitation, however, can contaminate 
brightness temperatures. Operational algorithms usually flag rain-contaminated 
radiances for rejection. 

Although the arctic regions have low humidities, absorption by water 
vapor can effect remotely sensed emissions. For the SSM/I only the 22 and 85.5 
GHz channels are significantly effected by water vapor absorption (Hollinger et 
al, 1987). Neither channel is currently used in ice studies, although the higher 
resolution offered by operational 85.5 GHz channels in future DMSP satellites 
will be tempting. 

The atmospheric oxygen distribution is well Known and constant. 
Absorption effects are included in seasonal climatological constants used in ice 


concentration algorithms. 


C. SEA ICE CONCENTRATION ALGORITHM 
1. Physical Model 


A simple radiative model expressing the radiation sensed at satellite 


height (Tp) is presented in Svendsen et al. (1983). Surface and atmospheric 


effects are included, plus a minor free space contribution. The model is composed 
of four terms (see Figure 3): radiation from the surface as seen through the 
atmosphere (A); upwelled radiation from the atmosphere (B); downwelled 
atmospheric radiation reflected at the surface (C); and free space radiation 


reflected by the surface (D). 


TRB = Esfel sfc(1-Ta) + OTaty + (1-Ec6-)6T,Ta(1-T.) + (1-€cfo)T sp (2) 
(A) (B) (C) (D) 


Tq is the optical depth, Tsp is the free space temperature and Ta is the weighted 
average atmospheric temperature in the lower troposphere. Since the optical 
depth is small, the extinction exponential (e-t4) is approximated by 1-Tgq. 
Atmospheric absorption of free space radiation is neglected. 

From equation (1) the emitted brightness temperature (Te) is dependent 
only on the physical temperature of the emitter (Tsf-) and its emissivity (€sfc). By 
substituting Te/Tsfc for sfc in equation (2), Te can be computed if tg and 6Tg are 
known, and if Tsf¢ is inferred from 5T3. 

Three components contribute to Te: sea water, first year ice and 
multiyear ice. The sum of the fractional components of each (Cy, Cf and Cm, 


respectively) add up to 1: 
Ee Ceo Ore Orr (3) 


The emitted brightness temperature is the sum of the individual components, 


weighted by the concentrations: 
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Te = Este l ste = CmEmTice + Cet lice + Cuew2/2) (4) 


where water temperature is assumed to be 272° K and €f, €m and Ew are the 
emissivities of first year ice, multiyear ice and seawater, respectively. Tice in 


equation (4) is determined by assuming a linear relationship between the 


atmospheric surface temperature (T.,) and the water temperature: 


Nea Ol ce lea) 72) ©) 


where & is an observationally determined constant (= 0.4 from Svendsen et al., 
1983). 

Fractional concentrations can be computed from equations (1) - (5) 
using two radiometer channels. Sources of error include approximations of 
optical depth, ice temperature, atmospheric temperature weighting factor, 
atmospheric surface temperatures and assuming constant emissivities for first and 
multiyear ice. Svendsen et al. (1983) compared aircraft radiometric measurements 
with Nimbus 7 radiometer measurements using this technique with in-situ 
temperature observations and determined accuracies to be + 3% of total ice 
concentration and + 10% for multiyear ice concentration with a spatial resolution 
of 90 km. 

2. Multichannel Algorithm 
a. Multispectral Analysis 
Comiso (1983) applied a multi-spectral technique that exploits the 


variation of emissivity with frequency for sea water and various ice types to 


I] 


distinguish between various first and multiyear ice and open water. Figure 4 
shows three distinct brightness temperature clusters that, when geo-located, 
correspond to first year ice, a first/multiyear ice mixture and open ocean 
conditions. The line that connects the first year cluster with the open ocean is 
interpreted as varying ice concentrations. Because the skin depth is frequency 
dependent, the grouping beneath the first year ice is theorized to be thin ice, 
where cold contamination from sea water is observed in the lower (deeper 
penetrating) frequency. The larger spread of emissivities and brightness 
temperatures in the 37V channel reflects the greater sensitivity of 37 GHz 
radiation to near-surface effects due to its smaller skin depth. 
b. Polarization/Gradient Ratios 

Cavelieri et al. (1984) used polarization ratios (PR) for a single 

frequency and spectral gradient ratios (GR) between two frequencies to estimate 


ice concentrations and fractions. PR and GR are defined as follows: 


Tpyfy) - Tpn(f) 


PRS Tal) + Tay . 
_ Tpy(fy) - Tay{la) 
R= Teh) + pv) Mi 


where Tp, and Tp, are brightness temperatures for vertical and horizontal 
polarizations respectively and f, and f> are separate frequencies. Table 2 shows 
PRs and GRs for sea ice and water emissivities at the two SSM/I frequencies (37.0 
and 19.35 GHz) used to calculate ice concentrations and fractions. Note that PR 


is much greater over the open ocean than over sea ice. Thus PR is a good 
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water/ice discriminator and can be used to estimate ice concentrations. GR 
depends on ice type and is affected by season. Ice types are determined by the 
sign of GR: positive over the open ocean, negative over multiyear ice and near 


zero over first year ice (Parkinson and Cavalieri, 1989). 


TABLE 2: EMISSIVITIES AT 19.35 AND 37 GHZ 
First Year | First Year | Multiyear Calm 
Thin Ice Ice Ice Sea Water 


Frequency 
(GHz) 


Sf 


| 379 [096] 087 | 097] 095 | 069] 0.64 | 070] 040 
Polarizati 

Teel 0.02 0.05 0.30 
(37V - 37H) 

erie 

= prea er 
(37V - 19V) 


Cavalieri et al. (1984) describe how these ratios are used to 





calculate total and multiyear ice concentrations and the physical temperature of 
the sea surface. They conclude that the precision of the calculated ice 
concentrations ranges from 5 - 9 % and the precision of the calculated multiyear 


ice fractions vary from 13 - 25 %. A caveat on the ice fraction calculation is that 


lig 


first and multiyear ice can not be discriminated during the summer season 
because the general surface melt conditions mask the distinguishing properties of 
the two ice types. A non-iterative version of this method has been adopted by 
the NASA Sea Ice Algorithm Working Group and is used by the NSIDC (NSIDC, 
1990). 


D. DATA SOURCES 
1. NSIDC CD-ROMs 

NSIDC distributes SSM/I brightness temperatures on CD-ROM disks for 
both arctic and antarctic regions. A continuous record is currently available from 
9 July 1987 - 31 December 1989. Earlier SMMR data are available in the same 
format from 1978-1987. Subsequent SSM/I data are in production for CD-ROM 
distribution. The data include all SSM/I channels with 85 GHz data gridded at 
12.5 km resolution and 19, 22 and 37 GHz data gridded at 25 km resolution. 
Added features are land masks and coast outlines that provide geographic 
reference points on the imagery. The data are obtained from 
NOAA/NESDIS/SDSD. Fleet Numerical Oceanography Center produces the 
original TDR files from the raw satellite data as part of their operational data 
processing. The TDRs are transferred to NOAA/NESDIS/SDSD and then to 
NSIDC. At NSIDC the TDRs are examined for bad data values and a geo- 
location correction is applied. The antenna temperatures are then converted to 
brightness temperatures. 

Swath data from each channel are mapped onto the appropriate grid by 
converting the geodetic latitudes and longitudes at the center of each 
observation footprint into map grid coordinates. Brightness temperatures are 


summed in each grid cell over a 24 hour period and then divided by the number 
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of observations. The resulting data fields are daily averaged mosaics of each 
polar region. 

The data are mapped onto polar stereographic projections with the grid 
tangent to the earth at 70° latitude rather than at the poles. This increases the 
distortion at the poles by 3% but is offset by nearly true projections in the 
marginal ice zones. Figure 5 shows the area of coverage for the Northern 
Hemisphere. 

2. Ice Buoy Tracks 

The Arctic Buoy Program has been operating since 1979. Objectives of 
the program are to provide surface atmospheric pressure measurements over the 
arctic basin and to define the large scale field of sea ice motion (Thorndike and 
Colony, 1980). Buoys are parachuted from C-130 aircraft at designated 
locations. Each buoy transmits pressure and air temperature information to a 
satellite at one minute intervals. This information 1s only received when the 
satellite is within the radio field of view. The data are retransmitted to a ground 
station, decoded and filtered, and forwarded to the Polar Science Center at the 
University of Washington. Thorndike and Colony (1980) examined the 
locational record of one buoy and reported the positional accuracy, as measured 
through the Doppler shift of the buoy signal received by the satellite, to be within 
300 m for two thirds of the reported positions. Very few fixes were reported with 
errors greater than 10 km and these can be easily removed from the record by 


careful editing. 
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Figure 1 Schematic illustration of differences between (a) first year ice and 
(b) multiyear ice (after Comiso, 1983). 
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Figure 2 Emissivities of calm seawater, multiyear ice and first year ice as a 
function of frequency for horizontal and vertical polarizations 
(modified from Svendsen et al., 1983). 
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Figure 3 Illustration of the four terms in the radiation model. 


18 


Absorbtion 





(b) SMMR Brightness Temperatures (K) 
February 25, 1979 
Section 1 


Zia 


Ta (V8 GHz, V) 


7 
i 
) 
1 


Tg (37 GHz, V) 





Figure 4 Multispectral analysis of brightness temperatures (from Comiso, 
1983). 


19 





Figure 5 SSM/I north polar grid (from NSIDC, 1990). 


Ill. CASE STUDY: ANOMALOUS DIURNAL CURRENTS ON 
THE YERMAK PLATEAU 


A. BACKGROUND 
1. The Yermak Plateau 

The Yermak Plateau (Figure 6) is a hooked, finger-like extension of the 
Barents Sea continental slope, located off the northwestern coast of Svalbard. 
The plateau extends 400 km into the Arctic Ocean and covers roughly 20,000 
km. Plateau depths range from 1000 to 600 m. The shallowest depth is 425 m 
on a ridge striking north-northwest near the plateau's central eastern edge (Perry, 
1986). The plateau slopes down rapidly to 3500 m to the north and northwest. 
The eastern slopes gradually level out to a semi-enclosed 2000 m deep basin. 

2. Regional Tides 

Tides in the Arctic Ocean are predominantly semi-diurnal. Schwiderski 
(1986) developed a model for global tidal prediction based on a 1° by 1° grid 
with accuracies better than 10 cm. Figures 7 and 8 show the structure of the 
leading semi-diurnal (M2) and leading diurnal (K1) partial tides, respectively, for 
the Nordic Seas. Dashed lines depict the co-range lines measured in centimeters. 
The M2 amplitude varies from zero at the amphidromic points to 80 cm and 
ranges from 20-30 cm over the Yermak Plateau. The amplitude of the K1 tide 
varies only from zero to 10 cm, ranging from 2-5 cm over the Plateau. 

Gjevik and Straume (1989) developed a regional tidal model to simulate 
the M2 and Ky tides for the Arctic Ocean and Nordic Seas with a 50 by 50 km 


grid size. The M2 and K1 amplitudes are similar to Schwiderski's global model. 


ZA 


The regional model also shows clockwise K1 current rotation along the west 
Norwegian shelf up to Svalbard. These authors believe this indicates the 
presence of a coastal trapped wave modal structure. 
3. Observations 

Anomalous diurnal currents over the Yermak Plateau were first reported 
by Hunkins (1986) from FRAM III and IV current meter data. The meters drifted 
from deep (4000 m) water to shallow water over the plateau's northwestern flank. 
Hunkins resolved absolute currents into orthogonal components in the along 
slope (60°) and cross slope (150°) directions and found that strong diurnal 
currents exist along the plateau's slope. Figure 9 shows current velocities at 104 
m along with the bathymetric profile for the FRAM IV drift. During the deep 
water portion of the drift, the current record is noisy with some evidence of semi- 
diurnal tidal motion. Over the slope diurnal motion dominates the signal. Once 
on top of the plateau, the semi-diurnal signal returns. Kinetic energy spectra 
confirm the dominance of diurnal over semi-diurnal amplitudes (Figure 10). 

Tidal motion polarization differs for the two tidal periods. Figure 11 
shows a broad diurnal peak in the kinetic energy spectra near 150° and a 
minimum near 60°. Conversely, the semi-diurnal motions are nearly parallel to the 
slope with a maximum at 60° and a minimum at 150°. The semi-diurnal alignment 
agrees with the modeled M? tides for the Arctic Ocean (Kowlik, 1979 as cited by 
Hunkins, 1986). Hunkins believed the cross-shelf diurnal currents indicated 
strong topographic vorticity waves. 

Strong diurnal currents were also reported from MIZEX 84 data. 
Morison et al. (1987) noted semi-diurnal currents dominated over most of the 


current meter drift record with amplitudes of 10 cm/s. Diurnal motions, however, 


Ze 


predominated when the current meter drifted over the southwestern slope of the 
Yermak Plateau. The diurnal component of motion ranged up to 12 cm/s. 

The observations of Hunkins (1986) and Morison et al. (1987) were 
confirmed by Hoffman's (1990) analysis of ARCTEMIZ 86 drifter buoy tracks. 
He showed a strong correlation between the enhanced diurnal currents and 
bathymetry. While in deep water no significant peak in the energy spectrum at 
the semi-diurnal or diurnal period was found. However, once the buoys crossed 
the 3000 m contour into shallower water the energy spectra exhibited a diurnal 
peak which exceeded the semi-diurnal peak. Figure 12 is a buoy track that 
clearly shows diurnal clockwise looping over the plateau's northern slope. 
Hoffman also reported a periodic increase in diurnal clockwise motion that he 
associated with the fortnightly (spring-neap) beat in surface height inequality 
caused by the superposition of the two principal semi-diurnal tidal components. 

4. Analysis 

Hunkins (1986) showed that topographically trapped vorticity waves 
Can exist at the diurnal frequency over the Yermak Plateau. He ignored the 
details of forcing and used a simple free wave analytic model. Chapman (1989) 
explored the actual forcing of the waves by the diurnal tide. He showed that a 
rectilinear oscillating barotropic current is significantly amplified by Hunkin’s 
idealized representation of Yermak Plateau topography. 

Recent work by Prazuck (1991) indicated high kinetic energy cells are 
possible on the Plateau. His linear barotropic numerical model created such cells 
in close agreement with observed diurnal currents. He showed through an 


analytic model that the highest energy cell could be due to the backscattered 
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energy of a coastally trapped wave that is unable to propagate beyond the tip of 
the plateau. 

Prazuck concluded that the enhanced rotating diurnal currents may 
influence ice floe motion. The high energy cells would increase the grinding 
together of the floes and alter the break-up and freeze-up patterns in the region. 
Several phenomena might be related to this interaction: 

¢ Strong turbulent fluxes associated with the diurnal currents 
such as those found during CEAREX 89 by Padman and Dillon 
GE IU. 
¢ Enhanced ambient noise in the region of grinding. Cousins 
(1991) correlated strong variations in ambient noise levels with 
12 hour tidal/inertial oscillations east of the Yermak Plateau. 
¢ Altered ice morphology due to the convergence/divergence 
of ice floes. 

5. Effects on Sea Ice Characteristics 

The interaction of the enhanced diurnal currents on the sea ice should 
be a seasonal phenomena. The high ice concentrations of the winter pack ice 
should significantly dampen floe motion, reducing the possibility of altering ice 
conditions. Lower ice concentrations in summer on the other hand should permit 
larger amplitude motions and increase the chances of remotely observable 
changes of sea ice characteristics. 

Ice floe motion effects both ice concentration and morphology. The 
enhanced currents can induce local divergence in the ice field and raise or lower 
ice concentrations on a scale observable from space. Divergence and grinding 


together of floes can increase fracturing. Subsequent refreezing would increase 
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the fraction of first year ice in a predominantly multiyear ice field without 


significantly affecting the total ice concentration. 


B. METHODOLOGY 

Data from the NSIDC SSM/I CD-ROM data set were used to test the 
hypothesis that changes in ice concentration/morphology could be detected by 
Satellite remote sensing. SSM/I data were used because, although visible and IR 
measurements provide higher resolution data, cloud contamination impede most 
attempts to examine a time Series extending beyond a few days. 

Sea ice concentrations were computed using the algorithm adopted by the 
NASA Sea Ice Algorithm Working Group (NSIDC, 1990). Gradient ratios were 
also calculated by modifying the ice concentration program to return gradient 
ratios instead of ice concentrations. The 25 km grid spacing used by the NSIDC 
data set was retained. Appendix A contains the software code listing. 

Figure 13 shows how the Yermak Plateau region was divided into test and 
control areas. A strip of 5 pixels (125 by 25 km), extending over the 
northwestern flank of the plateau, was selected as the test area. Control areas of 
the same size were located 50 km to the northwest and 125 km to the southeast. 
These strips will be referred to as ‘Slope’, 'North' and ‘South’, respectively. A 
perpendicular strip 75 km upstream of the general ice flow, northeast of the 
plateau tip, was used to establish the gradient ratio of ice flowing into the test 
area. 

A 168 day time series of daily ice concentrations, extending from 1 May to 15 
October 1988, was studied. One four day and several one to two day data gaps 
in the record were filled by linear interpolation of surrounding values. As seen in 


Figure 14, ice concentrations varied from 75-100% during the period for the 
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northern and test strips, and dipped to as low as 35% for the southern strip. At 
other times of the year concentrations remained near 100%. For the gradient ratio 
analysis the time series extended from 13 January to 30 April 1988. This period 


avoided brightness temperatures contaminated by melt ponds. 


C. RESULTS 
a. Ice Concentration 

The ice concentration algorithm was first checked for spatial 
consistency. Figure 15 shows concentration values for each pixel within the 
Slope strip during the period of interest. Table 3 contains a correlation matrix 
showing that the values are spatially consistent with 0.74 as the lowest 
correlation between the two endpoint pixels. 

Differencing between the North - Slope and South - Slope 
concentration values (Figure 16) showed no significant deviations between ice 
concentrations over the slope and control areas. The North - Slope mean 
difference is 1.04%, well inside the 10% accuracy of the concentration algorithm. 
The Slope - South mean difference is larger at 8.56%, but as seen from the 
concentration record this is due to the generally lower concentrations in the 
South strip later in the summer. 

The record was next examined for evidence of the fortnightly beat of 
current intensity reported by Hoffman (1986). If the currents influence ice 
concentrations, the largest signal should be present when the currents are 
Strongest. An autocorrelation of Slope concentration values showed no 
Significant correlations over a fortnightly period (Figure 17). Spectral densities 
were also computed for the mean Slope, North and South concentrations, as well 


as the North - Slope and South - Slope differences (Figures 18 - 19). In none of 
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these cases was a spectral peak located at or near the 13.66 day lunar fortnightly 


cycle (0.07 cycles per day). 


TABLE 3: CORRELATION MATRIX FOR YERMAK 
SLOPE PIXELS 





b. Gradient Ratio 
The gradient ratio is indicative of ice age. As shown in Table 2, the 
gradient ratio should be near zero for first year ice, less than zero for multiyear ice 
and greater than zero for sea water. The Siberian Branch of the Polar Drift Stream 
advects mostly multiyear ice through this region (Polarstern Scientific Party, 
1988). If interactions between the diurnal currents and the ice are sufficiently 


strong to increase ice fracturing over the Yermak Plateau slope, the refreezing of 


ge) 


the resultant small leads should increase the total fraction of first year ice. An 
increased fraction of more saline ice will drive the gradient ratio more positive 
compared with the control areas to north, south and upstream of the flow. 
Gradient ratios for the North, Slope and South regions are shown in 
Figure 20. Mean values are -0.033, -0.031 and -0.034, respectively. All gradient 
ratios are negative as expected from the predominantly multiyear ice. Although 
the gradient ratios are at least 26% higher over the slope than over the control 
areas during days 10 - 42, the difference between upstream and slope values is 
negligible for the same period (Figure 21). This means that no change measurable 


by the SSM/I occurred as the ice flowed over the slope region. 


D. CONCLUSIONS 
Anomalous diurnal currents over the Yermak Plateau slopes have been 
reported and analyzed by Hunkins (1986), Morison et al. (1987), Hoffman (1990) 
and Prazuck (1991). These currents have the potential for affecting the sea ice 
cover in this region in several ways. Divergence/convergence patterns could 
alter local ice concentrations. Stresses induced in the ice could lead to increased 
grinding and fracturing of floes, changing the morphology of the ice cover. If 
such changes could be noted through remote sensing, existing data bases could 
be explored for similar phenomena in other ice-covered regions. 
Ice concentration/morphology changes are not observed by the SSM/I. Two 
explanations could account for this lack of confirmation: 
1) The currents have little or no effect on the ice cover. 
2) SSM/I limitations obscure features in the ice caused by the currents. 
Factors supporting the second hypothesis include poor SSM/I spatial resolution 


and contamination by melt ponds or snow cover. Passive microwave emissions 
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observed by the SSM/I are too imprecise to verify the anomalous current effects 


on the sea ice cover. 


Zo 





Figure 6 Bathymetric chart of the Yermak Plateau (from Hoffman, 1990). 
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Figure 9 Absolute current velocities along the FRAM IV track at 104 m (from 
Hunkins, 1986). 
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Figure 10 —‘ Kinetic energy spectra for the cross slope velocity component from 
FRAM IV at 303 m (from Hunkins, 1986). 
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Figure 11 Peak amplitudes of kinetic energy for the K, and Mo? frequencies as a 
function of current axis orientation for FRAM IV at 303 m (from 
Hunkins, 1986). 
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Figure 12 ein looping of ARCTEMIZ 86 drifter buoys (from Hoffman, . 
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Hatched boxes indicate 


Yermak Plateau test and control areas. 


individual pixels. 


Figure 13 


Depths are in meters. 
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Figure 14 Total ice concentrations, 1 May - 15 October 1988: (a) North 
control, (b) Slope and (c) South control. 
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Figure 15 Yermak slope ice concentrations, 1 May - 15 October 1988. 
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Figure 16 Percent concentration differences (1 May - 15 October 1988): (a) 
North control minus Slope, (b) South control minus Slope. 
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Figure 17 Autocorrelation of slope region ice concentrations. 
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Ice concentration spectra (1 May - 15 October 1988): (a) North 
Control, (b) Slope, and (c) South control. 
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Figure 19 —_Ice concentration spectra (1 May - 15 October 1988): (a) North - 
Slope difference and (b) Slope - South difference. 
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Figure 20 Gradient ratios for Slope and Control areas, 13 January - 30 April 
1988. 
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Figure 21 Slope - Upstream gradient ratio difference, 1 February - 30 April 
1988. 
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IV. CASE STUDY: ICE MOTION ANALYSIS 


A. BACKGROUND 
1. General Ice Flow 

Sea ice motion affects a variety of important arctic processes. Carsey 
and Holt (1987) note that ice movement affects latent heat advection, oceanic 
surface stress, open water production through divergence and shear, is a passive 
tracer of ocean currents, and affects icebreaker navigation and ice structure. 
Within the central ice pack, relative ice motion creates both open water, through 
which significant heat is lost, and ridges that are operational hazards. At the ice 
margins ice growth associated with divergence affects brine production which, in 
turn, influences intermediate and deep water formation. 

Ice motion is primarily due to local surface wind stress, ocean currents, 
Coriolis and inertial forces, and effects of horizontal stresses transferred through 
the ice itself (Hibler, 1986). Colony and Thorndike (1984) examined a 90 year 
record of ship, research station and buoy tracks to create a quantitative estimate 
of the mean ice motion in the arctic basin (Figure 22). Two features dominate the 
circulation pattern—the anticyclonic Beaufort Gyre in the Canadian Basin and 
the Transpolar Drift Stream (TDS) extending across the Eurasian Basin to Fram 
Strait. 

The TDS is the largest and most concentrated meridional ice flow in the 
world. It is the primary outflow region for arctic sea ice. Variations in flow 
magnitude are theorized to effect the global climate (Vinje and Finnekasa, 1986). 


As seen in Figure 23, the TDS is composed of two branches. The Polar Branch 
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advects old multiyear ice across the pole and out through Fram Strait. The 
Siberian Branch contains thick, deformed multiyear ice floes that probably 
originated on the Siberian shelves. Closer to the Yermak Plateau, generally thin, 
first year ice formed locally or in the Barents and Kara Seas is found (Polarstern 
Scientific Party, 1988). The merging of these different ice types makes Fram Strait 
a rich area for SSM/I feature tracking. 

2. Feature Tracking 

Feature tracking estimates motion of a fluid by tracking features that are 
advected in the fluid. Zhang (1991) identifies two main steps necessary for 
successful feature tracking: 

1) Tie point selection - Identification in the first image of suitable 
features to be tracked; and 

2) Tie point registration - Matching the same tie points in a subsequent 
image. 

The underlying assumption is that the tracked features are advected by 
the fluid. Motion vectors are representative of the surrounding medium only 
when this assumption is valid. Possible disruptive influences include waves in the 
fluid, when fluctuations in fluid properties are tracked rather than the Lagrangian 
flow of the fluid, atmospheric contamination which can temporarily obscure 
features, and Eulerian effects such as local heating or cooling that provide feature 
sources and sinks. 

In the past sea ice motion has been inferred from the drift records of ice- 
locked ships, manned research stations and buoy tracks. The number and spatial 


density of these in-situ observations is limited. Satellites offer the advantage of 
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repeated sampling over large areas—a potentially rich source of information on 
the spatial and temporal structure of the sea ice motion field. 

Sea ice movement can be computed from displacements of features in 
sequential images. Subjective techniques, which require manual tracking of 
specific floes, are labor intensive and subject to repeatability problems due to 
differences in operator training. Automated techniques, such as the maximum 
cross correlation method described by Ninnis et al. (1986), are much more efficient 
but must be heavily filtered to remove spurious movement vectors. 

Three electromagnetic regimes are currently available for passive space- 
based observations of ice motion. In order of decreasing resolution, they are 
visible (0.3 - 0.7 um), infrared (3 - 12 im) and microwave (0.3 - 1.5 cm). The major 
limitations of high resolution visible and IR imagery are their inability to penetrate 
cloud cover and the light dependence of visual wavelengths. For example, 
Emery et al. (1991) investigated Fram Strait ice motion with IR imagery but were 
limited to six images over a 10 day period in April 1986 that contained 
exceptionally cloud-free conditions. On the other hand, low resolution passive 
microwave imagery can penetrate clouds and are independent of light conditions. 

Active microwave systems offer the advantages of cloud penetration 
and light independence with much higher resolutions than passive methods. For 
example, Synthetic Aperture Radar (SAR) data from SEASAT had a resolution of 
25 m (Carsey and Zwally, 1986). SAR information from the ERS-1 satellite will 
soon be available, however its application will be limited due to the limited 
sampling in time and space, especially in the subpolar regions (Zhang, 1991). 

Many researchers believe satellite passive microwave sensing 1s 


unsuitable for tracking ice motion (Zhang, 1991; Cavalieri, 1991 personal 
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communication) due to its inherent coarse spatial resolution. The primary use of 
this satellite data for arctic sea ice studies has been to examine the temporal and 
spatial variability of the ice cover. Svendsen et al. (1983) reported the motion of 
a patch of high multiyear ice concentration that correlated well with buoy 
movements, but the potential for feature tracking has not been pursued in the 
published literature (Cavalieri, 1991 personal communication). 

One can understand why passive microwave feature tracking has been 
overlooked. As mentioned earlier, brightness temperatures are sensitive to 
variations in surface conditions such as wet snow or melt ponds. The spatial 
resolution is too coarse to track individual floes. Judicious application of SSM/I 
data, however, can overcome these limitations. The basic requirements are large 
scale spatial gradients in ice type (as determined from gradient ratios) to overcome 
the coarse data resolution, and small temporal variability of surface conditions to 
limit melt pond or wet snow contamination. The first requirement is achieved by 
examining areas where old and new ice merge. The second requirement is met by 
rejecting data during the summer melt season. 

These two restrictions limit the wide-scale use of inferring ice motion 
through SSM/I feature tracks. Regions where this technique shows the most 
promise, however, are areas of high interest such as Fram Strait, through which 
occurs approximately 90% of the heat exchange and 75% of the mass exchange 


between the Arctic and World Oceans (Aagaard and Greisman, 1975). 


B. METHODOLOGY 
1. Satellite Images 
Multiyear ice fractions were calculated using the NASA gradient ratio 


technique described earlier. The same software that computed ice concentrations 
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and gradient ratios for the Yermak Plateau case study was used to create 
multiyear ice fraction images. Appendix A contains the FORTRAN code listing. 
The software merges the computed multiyear ice fractions with a land mask to 
produce a 304 by 448 pixel polar stereographic representation of the arctic 
region, extending from the pole south to 70° N. The pixel resolution is 25 km. 

Images were produced from the NSIDC Winter 88 and Winter 89 data 
sets. Images were created at 3 day intervals. If one assumes a typical 10 cm/s ice 
speed, a feature advected by the flow will move nearly 26 km over 72 hours. 
Since the grid resolution is 25 km, a 3 day sampling interval is the minimum 
needed to show expected movement of a feature within the data set. 

The imagery was displayed on a Macintosh II desktop computer using 
Imagic image processing software (Norris and Emery, 1989). Imagic was designed 
to display satellite remote sensing imagery on the Macintosh and is distributed by 
NSIDC with the CD-ROM data set. Images were color enhanced to accentuate 
multiyear ice features and looped to provide a continuous animation of feature 
movement. 

2. Feature Tracking 

The two main steps for determining ice motion from successive images 
are tie point selection and registration. These steps are achieved in this study 
through subjective imagery analysis following rules that could be incorporated 
into an objective scheme. Tie points are selected by identifying spatially coherent 
features through a visual inspection of the animated imagery. These tie points are 
registered in succeeding images by matching the same points in their new 
locations. This procedure is described in detail below: 


1) Create a series of multiyear ice fraction data sets. 
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2) Load the data sets into Imagic and generate ice fraction images. 
Decide on a color enhancement scheme that best shows multiyear ice features 
and apply this enhancement to all images. 

3) Using Imagic, observe the advection of the multiyear ice features. 

4) Select features that appear to be advected with the general ice flow 


according to the following criteria: 


¢ Size - Features must be small, no more than 100 km across (4 
pixels. Accurate tiepoint registration within larger features is too 


difficult. 


¢ Background Discrimination - Features must stand out from the 
background by at least 10 percentage points for most of their 
lifespan. This increases one's confidence that what is being 
observed 1s actually a feature advected in the flow rather than 


random variability of multiyear ice fractions. 


¢ Temporal Coherence - Features must be present in a least 6 
consecutive images (18 days). This further reduces the possibility 
of spurious tracks due to temporary brightness temperature 
fluctuations in the ice cover. 
5) Record the pixel position of an appropriate tie point within the 
feature during its lifespan using the following guidelines: 
¢ If one pixel value is 3 percentage points or more higher than all 


other values within the feature, select it. 
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¢ If 2 or more pixel values are within 3 percentage points, select 
the pixel which best corresponds to the previous relative position 
of the tie point within the feature. 

¢ If the relative maximum within the feature is ambiguous, 
choose the pixel which most conserves ice motion. 

6) Convert the pixel location within the polar stereographic grid to 
geodetic latitude and longitude using software supplied by NSIDC. A version 
modified for the Macintosh ts listed in Appendix A. 

7) Determine ice velocities by measuring the change of position with 
time. FORTRAN code that converts geodetic positions to u and v components of 
velocity 1s also contained in Appendix A. 

An example of the feature tracking technique is shown in Figures 24 
through 26. This particular feature was tracked from 1 February through 22 April 
1988 as it moved southward down the east coast of Greenland. Figure 24 
contains six consecutive images of multiyear ice fraction contour plots from 3 - 17 
March. Each plot covers the same 500 by 250 km region centered at 76.1°N, 
13.5°W (note the Greenland coast in the lower left corner of the plots). An area 
of relatively high multiyear ice fractions can be seen migrating from the upper left 
to the lower right corner of each plot. This is the trackable feature. Figure 25 
shows pixel values for this same region expressed as multiyear ice fractions. The 
circled number is the feature tie point. Pixel locations were converted to geodetic 
coordinates and the u and v velocity components are shown in Figure 26. 

3. Feature Selection 
The number of trackable features depends on the region of interest. For 


example, the generally homogeneous conditions in the high arctic (above 85°N) 
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Create a poor environment for feature tracking. No suitable features were 
observed in the 1988 data set. Only two features in the high arctic were 
observed during the 1989 winter season. On the other hand, Fram Strait is a 
potentially rich area due to the merging of the Siberian and Polar branches of the 
TDS. 

To assess the validity of this feature tracking technique, the trajectories 
of multiyear ice features were compared with true ice motion using the arctic data 
buoy program as ground truth for the ice velocities. The criterion to colocate 
trackable features with a drifting buoy was challenging. Suitable features had to 
be identified and their positions compared against buoy tracks as reported in 
Colony and Rigor (1990 and 1991). Three feature-buoy matchups were selected 
from the winter 88 data set, each representing a different region of the arctic. 
Named after the closest geographic reference point, the features are Severnaya 
Zemyla (SZ), Franz Josef Land (FJL), and East Siberian Sea (ESS). A feature from 
the Winter 89 data set (CA) was chosen to represent motion in the central arctic 
basin. Although no satisfactory matchup between a trackable feature and a 
drifting buoy was found in the Fram Strait region, a buoy (#7001) that passed 
down the Greenland coast two weeks ahead of a feature (FRAM) was used for a 
relative comparison of the ice trajectory. 

The tracked features and ground truth buoys were not coincident in 
space; the average distance between the selected features and buoys varied from 
76 to 557 km. The closest matchup (76 km) occurred between FJL and buoy 
#1790. The ESS feature was 420 km from buoy #1896. SZ was 320 km from the 
buoy #7048. CA averaged 497 km from the closest buoy (#3288). Table 4 


summarizes the characteristics of the feature-buoy tracks. 
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TABLE 4: MYI FEATURE - BUOY DATA 


(Total Days 


845N. |31MAR89|_ 86.6N. 
3288 1JAN89 | 443.55 (90) 125.2E 
79.9N, |22APR88] 68.1N. 


Geolocation Key: - Severnaya Zemlya Note: Features are named 
a L - Franz Josef Land according to nearby 
ESS - East Siberia Sea geographic reference 
CA - Central Arctic points. 
FRAM -Fram Strait 





4. Errors 

Before ice motion vectors between the buoy - feature pairs are 
compared, the maximum errors inherent to this technique are described. Three 
types of positional errors are possible: imperfect tiepoint registration, SSM/I data 
geo-location, and ARGOS buoy positioning. In addition, ice velocity gradients 
between the feature and buoy locations need to be considered. These errors, 
expressed as an ice speed error over a 3-day period, are examined in detail below. 

Imperfect tiepoint registration within the data grid is the largest source 
of error in this feature tracking technique. For small, well-defined features the 


maximum registration error is estimated to be one pixel for each observation. The 
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actual error should be uniformly distributed within the pixel. The maximum error 
is equivalent to a 50 km or 19.3 cm/s velocity error over a 3 day period. The geo- 
location error, as reported by NSIDC (1990), is 8 km for each scan. Since each 
pixel represents a daily averaged mean of multiple scans, the geo-location error 
Should be somewhat smaller. For an error of 5 km, the magnitude of ice motion 
could vary by 3.8 cm/s over 3 days. Thus the maximum estimated position error is 
30 km. If this error is assumed to be uniformly distributed, then the standard 
deviation of the error is 19.0 km. 

The estimated velocity error due to the maximum position error is 23.] 
cm/s over 3 days. This is over twice the typical ice velocity of 10 cm/s. If instead 
the standard deviation of the position error is considered, the velocity error 
reduces to 14.7 cm/s. These velocity errors decrease if the motion is averaged 
over longer time periods. For example, a 12 day average would result in a 
maximum error of 5.8 cm/s. 

The ARGOS system currently used to track buoy locations can 
determine position accuracy within a few hundred meters (Thorndike, 1986). 
Assuming a maximum error of 500 m for each ARGOS reported position, the 
resulting velocity error over 3 days would be 0.4 cm/s. 

When comparing SSM/I-derived ice trajectories with buoy locations, the 
velocity gradient between the two locations must be considered. These velocity 
gradients are the result of the ice pack response to external forces such as wind 
and ocean currents. Thorndike (1986) compared the velocity components of ice 
camps and buoys separated by distances ranging from 1 to 1000 km. He 
concluded that on scales of this magnitude ice behaves as a homogeneous, 


isotropic velocity field. Thorndike examined the correlations of buoy and ice 


ws 
A" 


camp velocities as a function of distance. For a 100 km camp-buoy separation, he 
found the correlation to be 0.9. At a 400 km separation, the correlation drops to 
approximately 0.6. 

Using these correlations, the variance of the differences in ice velocity 
between two separate points within the ice pack can be estimated. The variance 


of the velocity difference between the points (62 er) is given by 


O*err = E[V(x+) - V(x)]4 (8) 


where E is the expected value operator and V(x) and V(x+s) are the ice velocities 
separated by some distance s. Assuming the variance of ice velocity is the same 


at both points, (8) reduces to 


O-err = 207;(1-) (9) 


where 64; is the variance of the ice velocity and p is the correlation coefficient. 
Given a standard deviation of ice velocity of 10 cm/s (Thorndike, 1986) and 
assuming the u and v components of motion are independent of each other (the 
total variance then equals the sum of the variance of each component), O7¢,; is 40 
cm2/s2 and 10 cm2/s2 for correlations of 0.6 and 0.9, respectively. The resulting 
velocity errors are 6.3 cm/s for a 400 km separation and 3.2 cm/s for a 100 km 
separation. These velocity errors equate to 16 km and 8 km positional errors, 
respectively. 

Table 5 summarizes the maximum errors expected for each error source. 


The first three position errors are per observation. The corresponding velocity 
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error results from combining two position errors for each velocity calculation. 
Velocity errors are calculated over 3 and 12 day day intervals to show how the 


maximum estimated error decreases for longer mean flow calculations. 


TABLE 5: SUMMARY OF MAXIMUM ERRORS 


7 3 - Day 12 - Day 
Error Source a Velocity Error | Velocity Error 


(cm/s) (cm/s) 


Tie Point 95 
Registration 


SSM/I 
Geolocation 


vi 

Position 

Velocity 

Gradient 8 (3 Day) ceo) 

(100 Km) 

Velocit 

ened 16 (3 Day) 6.3 

(400 Km) 





This error analysis shows that 3 day ice motion estimates are unreliable 
indicators of true ice motion because the velocity error can exceed the typical 


true ice velocities of 10 cm/s. Longer term averages, however, reduce the impact 


>) 


of the SSM/I position errors and have the potential to represent accurately 
general ice flow. 

An expected correlation coefficient can be calculated from the maximum 
expected errors. This correlation can be used to validate the hypothesis that the 
feature tracks are representative of true ice motion. If the observed correlations 
are less than the expected value, then either the feature tracks do not represent 
true ice motion or the errors were underestimated. If the observed correlations are 
larger than the expected value, then the actual errors are smaller than expected. 


The correlation coefficient between two quantities, x and y, is given by 


Ex 
[Exx Eyy]!/2 


(10) 


where both x and y are assumed to have zero means. Let x and y be defined as 


the ice motion at the feature and buoy, respectively, as observed from the data: 


X=vt+eft+e€d (11) 
y=v+€b (12) 


where v is the true ice velocity, ef is the error in measuring the velocity of the 
feature, €d is the error due to the velocity difference between the feature and the 
buoy and €b is the error in measuring the velocity of the buoy. The variances due 


to these errors are defined as 


o2f = Eef2 + Eed2 (13) 
62h = Eeb2 (14) 
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62; = Ev2 (15) 


where of is the variance caused by the feature position and velocity difference 
errors, O~b is the variance due to the buoy position errors and 62; is the variance 


of true ice velocity. The expected values in equation (10) become 


Exy = 07} (16) 
Exx = 02) + o2f Gly) 
Eyy =07i+ 025 (18) 


Thus the expected correlation coefficient can be defined as the variance 
of true ice motion divided by the square root of the variances of buoy motion and 
feature motion multiplied together: 

G2: 


eee 
P [(o2, + 02, (02, + o2p 12 we 


As before, o?. for each component of motion is 50 cm2/s2. Sine tne oacay. 


sampling interval from Table 5, o- 4 is negligible at 0.16 cm2/s2 and o-- is the sum 


of the variances due to positional error (266.8 cm2/s2) and the velocity gradient 
between the buoy and the feature (39.7 cm2/s2 for 400 km and 10.2 cm2/s2 for 
100 km separations). Using these values in equation (19), the expected 
correlation for a feature and buoy 100 km apart is 0.391. For a separation of 400 


km, the correlation decreases to 0.374. 


These numbers indicate that if the correlation between 3 day mean 
velocities determined from buoys and multiyear ice features exceeds 
approximately 0.4, the feature is moving with the ice and the errors associated 
with the tracking technique were less than expected. If the same calculations are 
performed for 12 day mean ice velocities, correlations rise to 0.866 and 0.783 for 


100 and 400 km separations, respectively. 


C. RESULTS 
Il. Ice Trajectories 

Trajectories of the four feature - buoy pairs from the different arctic 
regions, each plotted at 3 day intervals, are shown in Figures 27 through 31. 
Because the feature tracks are limited to discrete 25 km jumps in position, periods 
of slow, looping or retrograde motion confuse the trajectories. In these cases, the 
positions have been slightly offset to show movement more clearly. 

The correspondence between 3 day positions of the buoys and features 
appears somewhat erratic. This is probably due to the resolution limitations 
inherent in the SSM/I data grid. Taken over longer periods, however, the general 
flow corresponds quite well to the buoy positions. This is best shown in the Fram 
Strait trajectories (Figure 27). Both the buoy and FRAM tracks are roughly 
parallel as they traveled southward down the east Greenland coast, the buoy 
preceding FRAM by about two weeks. Of special note is the sudden surge 
eastward at 78°N for the feature and at 69°N for the buoy. Although separated 
by 9° of latitude, this surge occurred at the same time (20-23 Feb). This 
corresponds to a period of strong geostrophic wind flow to the northeast due to a 
984 mb low centered over the northern coast of Greenland (Colony and Rigor, 


1990). 
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2. Velocity Comparisons 

For each feature - buoy pair, the u and v components of motion over 
each 3 day period were calculated. Buoy velocities were computed directly from 
their reported positions. With the exception of ESS, all multiyear feature 
velocities were smoothed over two positions to reduce velocity spikes resulting 
from the gridded data set. The ESS feature, tracked for 24 days, has only one 
local velocity maximum in each direction so the smoothing was not applied. 

Correlation coefficients between the buoys and both smoothed and 
unsmoothed feature velocities are presented in Table 6. Applying at distribution 
test to the smoothed velocities, all but one correlation (Franz Josef Land u 
component) exceed the 95% confidence level that the feature and buoy 
velocities were correlated. The generally high correlations are particularly 
Surprising since the maximum error analysis indicates that correlations only as 
high as 0.4 should be expected for 3 day mean ice velocities. These higher than 
expected correlations imply that the estimated errors used in the error analysis 
were too large. 

Specific velocity comparisons are plotted in Figures 32 through 35 and 
are discussed below. 

a. Franz Josef Land - Buoy 1790 (Figure 32) 

This is the closest spatial match-up of any feature - buoy pair 
investigated. The separation distance averaged 76 km over a 42 day period. The 
buoy trajectory shows complex motion with two separate periods of retrograde 
east - west movements. The SSM/I grid had trouble resolving these complicated 
movements, providing the lowest correlation of any velocity component: -0.069 


for the unsmoothed and 0.322 for the smoothed u component. When velocities 
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were averaged over 12 days, the unsmoothed u component correlation improved 
to 0.666. Three day v correlations were 0.458 and 0.678 for unsmoothed and 
smoothed velocities, respectively. The higher correlations in the v direction are 


probably due to the slow, southward trend of the drift. 


TABLE 6: MYI FEATURE - BUOY CORRELATIONS 


Feature Correlations (3 day velocities) 
Kael 


724 845 490 





b. East Siberian Sea - Buoy 1896 (Figure 33) 

This pair was tracked over the shortest interval (24 days). Although 
separated by 320 km, it also had the highest correlations. The correlations for the 
unsmoothed u and v components were 0.791 and 0.724, respectively. The 
trajectory characteristics responsible for these high correlations were a steady 
northwestern drift interrupted by a strong southwesterly surge (16-19 April) large 
enough to be resolved within the SSM/I grid. This surge corresponded to a 
period of strong geostrophic wind flow to the southwest from a low in the 


Chukchi Sea region (Colony and Rigor, 1990). 
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c. Severnaya Zemyla - Buoy 7048 (Figure 34) 

The feature and buoy were tracked for 75 days. The average 
separation distance was 420 km. Both trajectories showed a slow easterly drift 
culminating in a clockwise loop. Correlations for the smoothed u and v 
components were 0.685 and 0.455, respectively. The higher u correlation is due 
to the predominantly eastward movement of the feature and buoy. 

d. Central Arctic - Buoy 3288 (Figure 35) 

This pair was taken from the 1989 data set to represent ice motion in 
the central arctic region. They were tracked for 90 days at an average separation 
of 491 km. The smoothed u correlation was low at 0.311. The v component was 
much higher at 0.582. The low u correlation can be explained by the steady 
northerly movement of the feature. East-west velocities were generally small and 
not well resolved in the SSM/I data grid. Correlations for the 12 day mean 
velocities jumped to 0.783 and 0.776 for u and v components, respectively. The 
longer sampling period allowed time for the slow motion of the ice field in this 
region to be resolved within the SSM/I grid. 

3. Position Error 
The maximum estimated feature position error can be compared against 


the position errors calculated from the data. Solving equation (19) for oF: 


l 
o7¢ = Oo, (2 - ] (20) 


o2; is assumed to be the same as the variance of buoy velocities and o2b is 


negligible compared with o7ice. The average variance of each buoy velocity 
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component is 30.08 cm2/s2 and the average correlation (p) between the 
smoothed feature and buoy velocities 1s 0.545. 

Using equation (20) and the average variances and correlations 
observed from the feature track data, o2f is 71.19 cm2/s2. The ice velocity 
difference between the buoy and feature must be taken into account because 62f 
is the sum of the variances due to position error and the velocity difference. 
Substituting the average standard deviation of the buoy velocities into equation 
(9) and again assuming that u and v components of motion are independent of 
each other, the variance due to the ice velocity gradient is 12.32 cm2/s2 and 3.08 
cm/s? for 400 km and 100 km feature-buoy separations, respectively. Taking 
into account the average feature-buoy separation was 326 km, nearly 17% of the 
observed variance can be explained by the ice velocity gradient between the 
feature and the buoy. The remaining variance of approximately 58 cm2/s2 
equates to a 10.7 cm/s total velocity error (u and v velocity variances added 
together). 

This velocity error is less than half of the maximum expected value of 
23.1 cm/s and 72% of the standard deviation of the expected error. A 10.7 cm/s 
velocity error over a three day period implies close to a 15 km position error in 
each of the two observations used to calculate the velocity. This is less than both 
the maximum expected postion error (30 km) and the standard deviation of this 
error (19 km). Possible explanations for the difference between expected and 
actual values include: 1) the tiepoint registration and geo-location errors are 
smaller than expected because the gridded brightness temperatures in the dataset 
are obtained by averaging information gathered from several satellite passes each 


day; 2) the actual ice velocity gradient between the features and buoys is smaller 
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than Thorndike’s (1986) estimates; and 3) the standard deviation of the ice 
velocities is smaller than the assumed value of 10 cm/s. The latter explanation 
probably had the greatest influence because standard deviations of the buoy 
velocities only ranged from 3 to 9 cm/s. These smaller than expected standard 
deviations could result from the the use of three day ice velocity averages which 


reduce the effect of short term velocity surges. 


D. DISCUSSION 

There is good reason to doubt the practicality of inferring sea ice motion by 
tracking multiyear ice concentrations from space-based passive microwave 
radiometers. Data resolution is poor. At 37 and 19.35 GHz (the frequencies 
currently used to calculate ice fractions), resolutions are 33 and 56 km, 
respectively. The 25 km grid resolution of the SSM/I dataset is only achieved by 
averaging several observations over the same spot each day. Not only does the 
low resolution increase the position error, it also greatly reduces the possibility of 
tracking individual floes. The radiometer must instead observe broad regions of 
anomalous radiances to create trackable features. The precision of the multiyear 
ice fraction algorithm is only estimated to be in the range of 13 - 25 %. Moreover, 
a variation of 0.04 cm in cloud liquid water content will contribute as much as 
30% to the uncertainty of computed multiyear ice fractions (Cavalieri et al., 
1984). 

In spite of these handicaps, ice drift inferred through tracking multiyear ice 
features from SSM/I radiometer data corresponded quite well with ice motion as 
measured by in-situ buoys. Ice feature trajectories were similar to nearby buoy 
drift patterns. In two cases (ESS and FRAM), when high atmospheric pressure 


gradients caused anomalously strong wind stress on the ice pack, both the feature 


and buoy exhibited similar trajectories in accordance with the forcing. 
Comparisons between the orthogonal velocity components of the feature and 
buoy motion averaged over 3 days also showed surprisingly high correlations. 
Instances of high correlations in one component direction but low correlations in 
the other can be explained by slow or oscillatory movement in the direction of 
the velocity component exhibiting the low correlation and are a function of the 
low spatial resolution of the grid. 

The good results achieved above were the result of being extremely selective 
during the tiepoint selection process. First, all tracked features were of high 
multiyear ice fractions surrounded by a background of relatively low fractions. 
Multiyear ice is thicker than first year ice and therefore less subject to 
deformational and destructive processes that might obscure the characteristic 
radiometric signal. Second, the features were small—no more than 100 km across. 
This simplified the task of registering the same tiepoint in successive images. The 
variation of multiyear ice fraction percentages within a larger feature could 
preclude accurate tiepoint registration. Limiting the feature spatial extent also 
reduced the possibility of confusing an atmospheric disturbance with an ice 
feature. Third, imposing a temporal coherence requirement of 18 days ensured 
that the tracked object was not a result of random brightness temperature 
fluctuations in the ice cover. Finally, limiting the record to the winter season 
reduced the possibility of surface melt conditions contaminating radiances from 
the sea ice. Another benefit of the winter season is that moisture-laden clouds 
with their contaminating effects are generally restricted to the open ocean and 


the marginal ice zone areas. 
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The number of trackable features was limited. The abovementioned 
limitations restricted the number of suitable features to only a dozen or so during 
the winter season. Relaxing these restrictions should produce more trackable 
features (and possibly introduce spurious features) and is an area suitable for 
further research. Development of an objective scheme to identify and track 
features should also increase the number of available datapoints. 

Despite the above restrictions, this limited ability to observe and track areas 
of anomalous multiyear ice fractions can be useful for arctic studies. Ice fluxes 
from the central arctic have been shown to be the largest exporter of fresh water 
from the Arctic Ocean and the largest freshwater input to the Greenland Sea 
(Aagaard and Carmack, 1989). Available databases, extending back to 1978, can 
be used to assess the intra and interannual variability of ice transport. 
Researchers can identify previous, and possibly forecast future, pulses of ice 
types to examine their influence on regional heat and salinity budgets. Ice camps 
are generally located on relatively smooth, thick floes of multiyear ice for 
logistical and safety reasons. Potential sites are currently scouted by aircraft. 
This space-based technique offers the possibility of remotely identifying sites well 
in advance of deployment. As satellite instrumentation improves and emanations 
from higher microwave frequencies are utilized, surface resolution will also 


improve resulting in increased reliability of short term ice velocity vectors. 


E. CONCLUSIONS 

When images of arctic multiyear ice fractions calculated from winter season 
SSM/I data are examined in sequence, distinct features are observed to move in 
the general direction of climatological ice flow. These features exist throughout 


the arctic basin but are mainly located in regions where first year and multiyear 
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ice merge together. Some features are short-lived, others can be tracked 
throughout the data record. 

Four of these multiyear ice features from the winters of 1988 and 1989 were 
compared against nearby arctic buoy trajectories. Three day average ice velocity 
vectors were found to correlate well with the buoy velocities; the average 
smoothed velocity correlation was 0.545. The maximum estimated velocity errors, 
mainly due to the low data resolution, are double a typical ice velocity of 10 cm/s 
making short-term ice motion calculations unreliable. Actual velocity errors 
derived from the feature-buoy correlations are half the maximum estimated error 
and less than three quarters of the standard deviation of this error. This difference 
is probably due to the smaller than estimated variances of ice motion because the 
ice motion is averaged over a three day period, smoothing out short-term ice 
velocity surges. 

The higher than expected correlations between multiyear ice feature 
velocities and arctic buoys are due to the stringent requirements established for 
feature selection. Relaxing these requirements could substantially increase the 
number of trackable features from the dozen or so noted in the winter of 1988 but 
also decrease the confidence that tracked features accurately represent ice 
motion. A similar increase in trackable features could probably be attained by 
objectively scanning the dataset for subtle features that were not recognized 


through a subjective analysis of the data. 
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Figure 22. Field of mean ice motion (from Colony and Thorndike, 1984). 
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Figure 23 Average sea ice drift showing the two components of the 
Transpolar Drift Stream ( from Polarstern Scientific Party, 1988). 
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Iigure 24a Contour plots showing advection of multiyear ice feature along the 
east Greenland coast. Images from 2 and 5 March, 1988. Plots are 


centered at 76.1°N, 13.5°W. The feature moves from the upper left 
to the lower right corner. 
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Figure 24b Same as Figure 24a except images are from 8 and 11 March, 1988. 
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‘igure 24c ~=Same as Figure 24a except images are from 14 and 17 March, 1988. 
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Pixel values for multiyear ice contour plot shown in Figure 24a. 
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value is tiepoint registration location. 


Figure 25a 
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‘igure 25b Pixel values for multiyear ice contour plot shown in Figure 
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Pixel values for multiyear ice contour plot shown in Figure 24c. 


Circled value is tiepoint registration location. 


Figure 25c 
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Figure 26 Ice velocities for Fram Strait feature (1 February - 22 April, 1988). 
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Ice trajectories: Fram Strait feature and Buoy 7001 (1 February - 22 


April, 1988). Circled points are initial positions. 
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Figure 28 Ice trajectories: Franz Josef Land feature and Buoy 1790 (1 
February. - 14 March, 1988). Circled points are initial positions. 
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East Siberia Sea feature and Buoy 1896 (1 - 25 


April, 1988). Circled points are initial positions. 
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Severnaya Zemyla feature and Buoy 7048 (1 


February - 28 April, 1988). Circled points are initial positions. 


Ice trajectories: 
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Ice trajectories: Central Arctic feature and Buoy 3288 (1 January - 


31 March, 1988). Circled points are initial positions. 
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Figure 32. ~=Franz Josef Land and Buoy 1790 velocity comparison. Buoy is 
solid line, feature is dashed line. 
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Figure 33. ~~ East Siberian Sea and Buoy 1896 velocity comparison. Buoy 1s 
solid line, feature is dashed line. 
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Figure 34. Severnaya Zemyla and Buoy 7048 velocity comparison. Buoy is 
solid line, feature is dashed line. 
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Figure 35. = Central Arctic and Buoy 3288 velocity comparison. Buoy is solid 
line, feature is dashed line. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


Arctic sea ice conditions are an important component of the global climate 
system. The sea ice cover affects the Earth's radiative balance by influencing the 
heat flux between the ocean and atmosphere and heat transport between the 
lower and upper latitudes. Passive microwave remote sensing of sea ice offers 
continuous, all-weather observations of these remote areas independent of solar 
illumination. This is particularly important in the arctic where much of the region 
is shrouded by clouds or covered in darkness. 

Accurate ice observations using passive microwave remote sensing 
techniques require knowledge of the radiative characteristics of sea ice. The 
radiative properties of sea ice are complicated due to its complex nature which 1s 
determined by formation environment, age and season. Other complicating 
factors include contamination from moisture-laden clouds and the presence of 
snow cover, melt ponds and ridging. In spite of these difficulties, emissivity 
differences between sea ice and water at various microwave frequencies can be 
exploited to estimate ice characteristics. Current algorithms use polarization and 
gradient ratio differences to estimate ice concentrations and type. 

Historical SSM/I data are readily available on CD-ROM media distributed by 
NSIDC. Two applications of this dataset to arctic sea ice studies were examined. 
Data from 1 February to 15 October 1988 were used to associate anomalous 
diurnal tidal currents over the Yermak Plateau with ice morphology. Data from | 
February to 28 April 1988 and 1 January to 31 March 1989 were used to explore 
the possibility of determining sea ice motion by tracking advection of multiyear 


ice features. Software supplied by NSIDC were modified to compute and display 


the requisite ice concentration, gradient ratio and multiyear ice fraction 
calculations and imagery on a Macintosh IT desktop computer. 

The hypothesized ice concentration/morphology changes induced by the 
anomalously strong diurnal tidal currents observed over the Yermak Plateau 
could not be verified by the SSM/I data. Either the currents do not affect the ice 
cover or SSM/I limitations obscure the ice signatures of the currents. The low 
spatial resolution and susceptibility to melt and snow contamination of passive 
microwave emissions prevent a definitive assessment of tidal current effects on ice 
cover in this region. 

Distinct multiyear ice features are observed to move in the general direction 
of the mean ice flow. Feature movement correlates well with ice motion as 
measured by arctic data buoys. The subjective nature of the feature tracking 
technique and its reliability on ‘features of opportunity’ prevent wide scale 
application throughout the arctic region. Positional errors resulting from low data 
resolution could make short term ice motion calculations unreliable. Long term 
averages, however, are reliable and can be used to augment our knowledge of 
general ice flow. For example, motion vectors can be used to estimate the 
interannual variability of ice transport through the Fram Strait and aid our 
understanding of heat and salinity budgets in the Greenland Sea. In spite of the 
limitations, tracking multiyear ice features can be a useful method of discerning 
continuous, all-weather ice motion from satellite data. 

Further attempts to verify anomalous diurnal currents over the Yermak 
Plateau using SSM/I data are not recommended. The data resolution is too low to 
observe small scale changes of sea ice concentration or morphology. Feature 


tracking, however, does show some promise. Conversion of the subjective 


86 


technique described here to an objective scheme for computer analysis would 
greatly reduce the time and effort required for tie point selection and registration, 
and probably increase the number of tie points suitable for tracking. Exploration 
of the higher spatial resolutions afforded by 85.5 GHz data would greatly 


increase the reliability of ice motion vectors. 
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APPENDIX A 


This appendix contains FORTRAN source code listings of software used to 


process the dataset contained on the NSIDC CD-ROM disks. 


Program 


ssmlice 


locate 


UV Calc 


Subroutine 


usetie 


ice25 


Page 


89 


5 


102 


108 


Purpose 


Computes ice concentrations from 
SSM/I brightness temperature values on 
the NSIDC CD-ROM dataset. 


Transforms I, J coordinates of the NSIDC 
SSM/I polar grids to latitude and 
longitudes coordinates and vice versa. 


Calculates u and v components of 
motion from a list of user specified 
geodetic latitude and longitude 
positions. 


Purpose 


Computes coefficients for the ice 
algorithm in ssmiice using NASA- 
approved tiepoints. 


Computes ice concentrations from 
NSIDC CD-ROM data at 25 km grid 
resolution. 


program ssmuce 
eee cece Ce CCCCECTECCECCECCECEC CECE CCCECCECCE CCE CEC CEC ECCCCCEC CCC CCCCC 
G 
c SSMIICE - This program computes ice concentrations from SSMI TB 
grids residing on CD-ROM. 


Inputs 


SER DEFINED 
FAAFEEHEEFFEFEEEEFEFEEEEEEEEFEFEFEFEEEFEFEPTE+E FEF 


drive : Name of CD-ROM disk 


indate : Date in the format YYMMDD (eg., 870710) 
From this variable the CD-ROM is searched for the 
appropriate filename. The filename for the ice 
concentration algorithm is derived from this 
variable. 


ihem : Defines the appropriate hemisphere to apply this 
algorithm (1=North, 2=South) 


icepar : Defines the type of ice concentrations to derive: 
1=total, 2=multi-year, 3=first-year. 


INTERNAL 
FEFEEEEFEFHEEFEFEEEFEEFEEEEFEEEFFEEEEEEEF FAH E EE EF 


infile : The path, name, and descriptor string describing 
the data file from which the SSMI Tb data will be 
extracted from the CD-ROM. User must specify appropriate 
filepath. 
buff : Input buffer containing the SSMI Tb data read from 
the CD-ROM. Two lines of SSMI TB grid data will 
read with a single IO request. 


mask : Landmask used in computing the 25 km 
ice concentration grid. 


idiml : Number of parameters in an SSM/I Tb grid cell 
idim2 : Number of cells in the SSM/I Tb grid line 
idim3 : Number of grid lines in an SSM/I Tb grid 


grdtyp : Three character descriptor designating the grid type 
to access from the CD-ROM (N3B or S3B) 


LG ee > Si ek GP aL | GP TAK > Cam > Rl GP es Gam GA > Dae > Man, Ja, GP kee GP a ce, SP A P| Ale Sa GP a CB a el SAN Gh fed SSP a Ya > a pap tl P07 a ak dele > a PI Ga Bae” Ses! Go NN pe Ge ee 
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icetyp : Three character descriptor designating the extension 
for the resultant ice parameter image (N3C or S3C). 


€ : Array of tiepoint coefficients. Returned from USETIE 
and used in GETICE. Refer to NODS Handbook Appendix 
CC pp 23-24. 


Outputs 


INTERNAL 
FELEFEALEEFEFEEE FEF EE FEF EEF TEETH TE +++ +++ t44t+4¢+4+44444 


outfil : The name and descriptor string which reference 
the data file which will be used to store the 
ice concentration grid. This string 1s computed 
internally from the variables, indate and grdtyp. The 
file will be created in the default directory. 
Program will abort if file currently exists! 


The file dimensions will be 304 1-byte cell values 
with 448 grid lines for the grid type 'N3C' and 

316 1-byte cell values with 332 grid lines for the 
grid type 'S3C’. 


ice : Output buffer containing the computed ice 
concentrations for the 25 km total ice concentration 
grid. Totice is 316 bytes in length which corresponds 
to the length of an S3C raster line. 


Go) 0 8G) Ga Go, Cn C2 OG CeO Gy O30 1G Oo a OG «OC OOOO ©} 


c Vince Troisi 

c 5/3/1989 

C 

c 1/4/1989 - Modified for SSM/I CD-ROMs 

c 7/1/91 - Modified for Apple Macintosh and 25 km grid by W. A. Wright 

CCCCCCCECECCCCCCECCOCCCECECCEECCEECCCEE EC CEGECCenECecerececceececrcces 
integer idim1, 1dim2, idim3 


c 
c There are five channels in the N3B or S3B grid cells 
c 

parameter (idim1 = 5) 

parameter (maxlen = 316) 


integer lunit, Ounit, munit, olen,z 
integer ymd(6), ier, recln,ok 


character * 3 grdtyp, icetyp 
character * 9 drive 

character * 2 year, month, day 
character * 6 indate(95), ifile 
character * 3. mon(12) 
character *80 outfil,infile, mfile 
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integer * 2 buff(maxlen*idim]),abyte,bbyte 
integer *1 ice(maxlen) , mask(maxlen),atemp(2),btemp(2) 
real C12) 
equivalence (bbyte,btemp) 
equivalence (atemp,abyte) 
data mon 
far we, Pep WiAk APR MAY JUN JUL. AUG. SEP: 
+ ‘OCT, 'NOV','DEC;/ 


c Query user for search criteria 


c Query for name of CD-Rom drive 
print *, 'Name of CD-Rom?’ 
mead (. (a) )-drive 


c Query user for file dates 


5 print *, ‘Number of files to create? (95 max)’ 
print *, ‘Note: Read lasts approx. 10 min. per file’ 
read *,n 


c print *, ‘Enter file dates (yvymmdqd).' 
open (unit=29, file='Archs HardDisk:THESIS:daylist, 
+ status='old’) 
Ceol i=l 
read (29,'(a)') indate(1) 


10 continue 


pmnt*, these are the selected dates: 
do 20 1=1.n 
print *, indate(1) 

20 continue 


print *, ‘OK? (1-yes, 2-no, 3-stop) <ret>' 
read *, ok 


if (ok.eq.2) go to 5 
if (ok.eq.3) go to 6000 


c Query user for appropmiate polar region 


price, 

print *,'Select one of the polar regions: 
prnt *, 1. Norh: 

pant *, 2-South: 

prmt, enter lor 2) 
fead *. nem 
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c Query user for ice type to compute 
prnics 

print *,'‘Select type of ice characteristic to derive:' 

print *,' 1. Total Ice Concentration’ 

print *,’ 2. Multi-Year Ice Concentration’ 

pant ~, Enter Tons 

read *, icepar 


c Combine responses to form filenames and paths and to define 
c IO constants 


if (them.eq.1) then 


grdtyp = 'N3B' 

icetyp = ‘N3C' 

lrec = 3040 

idim2 = 304 

idim3 = 448 

npar =idim1*idim2 
else 

grdtyp = 'S3B' 

icetyp = 'S3C' 

lrec = 3160 

idim2 = 9316 

idim3 = 332 

npar = idim1*idim2 
endif 


do 5000 z=iin 
ifile = indate(z) 
read (ifile(3:4),'G2)'}mm 


infile=drive//':'//grdtyp//':19'//Afile(1:2)//:'//mon(mm)// 
tr '/Afile(1:6)//.'//grdtyp//; 1 


print *, ‘infile = ', infile 
c Open brightness temperature file on CD 


Init = 21 
open (unit=iunit,file=infile,status='old',access='sequential’, 
+  form='unformatted’,recl=lrec,iostat=ler) 
if (ier.ne.O) then 
print *,infile, ' not accessible! <ret>' 
pause 
Stop 
endif 
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c Construct output file name and open the file. 


outfil = ‘Archs HardDisk:THESIS:Conc Data:'//ifile(1:6)// 
+ '. /ficetyp 
pent: Output tile = ‘outfil 
ounit= 22 
open (unit=ounnt,file=outfil,status='new’ ,recl=idim2, 
+ 1ostat=ier,form='unformatted’, access='sequential’) 
if (ier.ne.O) then 
print *,'Error opening file ',outfil(25:35) 
stop 
endif 


c Open Landmask file 
mfile = drive//:TOOLS:'//grdtyp//MASK.DAT;1' 
print *, ‘mask file = ',mfile 
munit = 18 
open (unit=munit,file=mfile, access='sequential'’, 
+ form='unformatted',recl=idim2,status='old’) 

c Get ice algorithm coefficients derived from tiepoints 


call usetie (c,shem) 


C-------- Start 1ce concentration algorithm 
do 1000 Ip=1,1dim3 
c Swap bytes in buffer file for Apple Macintosh 


do 800 inb=1 ,npar 
read(iunit) bbyte 


atemp(2)=btemp(1) 
atemp(1)=btemp(2) 


buff(inb)=abyte 
S00 continue 
read(munit)(mask(inm),inm=1,idim2) 
c Average TBs into 25 km cells and determine ice concentrations 
call ice25 (buff,c,maxlen,idiim] ,idim2,1ce,mask,1cepar) 
c Wnite grid line of ice concentrations to file. 


write(ounit)(ice(iout),10ut=1 ,idim2) 


oie) 


1000. continue 
C-------- end ice algorithm 


close (unit=munit) 

close (unit=iunit) 

close (unit=ounit) 

close (unit=29) 
5000 continue 


6000 end 
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program locate 
G2GeCCe ECCCECCCECCCCCCCECEE EC CCE CCCE CECE CECE CCECCCCECECCCECOCCCCCCECCECC 


c 

c LOCATE - This program transforms I,J coordinates of an SSM/I grid cell 
@ to latitude and longitude coordinates. This program provides 

the inverse functions as well. LOCATE interfaces to the revised 
forms of the subroutines, MAPXY and MAPLL. 


User-defined Paramters: 


gtype : Integer supplied by the user to describe one of the three 
grid cell dimensions (12.5 km, 25.0 km, or 50.0 km). 


ihem_: Integer supplied by the user to describe one of the two 
polar regions (i=North , 2=South) 


itrans : Integer supplied by the user to describe the type of 
transformation LOCATE will perform (1=I.J-to-Lat,Lon; 
2=Lat,Lon-to-I,J) 


1,) : Integers supplied by the user when itrans = 1. These 
integers describe the position of a cell in an SSM/I grid. 


lat,lon : Reals supplied by the user when itrans = 2. These 
integers describe the latitude and longitude in an SSM/I 
grid which LOCATE will transform to an I,J grid cell position. 
Note: All latitudes and longitudes must be entered as 
positive numbers! 


Internal: 


x,y : Distance in kilometers from the origin of the grid 
fie] pole). 


alat, 
alon : Computed latitude and longitude returned from MAPXY. 


SGN — : Sign of the latitude (positive = north latitude, 
negative = south latitude) 


delta : Meridian offset for the SSM/I grids (O degrees for 
the South Polar grids; 45 degrees for the North Polar 
grids. 


kk —: Integer variable used for reorientation of the gnd. The 
grid is 'flipped' in the Y direction for transformations. 


SLAT ~ : Standard latitude for the SSM/I grids is 70 degrees. 


numy : Number of linesin an SSM/I grid. This attribute varies 
for each of the six grids. 


Cy air oe OO), (0) OO Ci OOO * OO) OC) OO a OG ta Oa Ca ONO OO) GO Oe Oee MeO), Gum (a. © 


cell.» Size of the SSM/iemd (.12,5.kim.25,0 Kae O10) ama) 


xydist : Distance from the origin of the grid in the cartesian plane. 
The x- y coordinates for ‘the edge of the lower left pixel 

is (3850.0, 5350.0) for the northern grids and 

(3950.0, 3950.0) for the southern grids. 


: Radius of the earth in kilometers. 

E __: Eccentricity of the Hughes ellipsoid 

E2 _ : Eccentricity squared 

Piel 

Written by V.J.Troisi - January, 1990 


Modified by W.A. Wright - July, 1991 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCECCCECCCECCCE CE CCC CCC CCC CC CCE CCC CCCEE CCC 


CIO) G25 (GG): OC). COG) OC) 6) C) OtG@ eG) CoC) 


real SLAT,E,RE,PI 

real alat,alon,x,y 

integer ihem 

real lat, lon 

real SGN, delta 

integer numy(2,3) 

real cells) xxvcisti ZZ) 

data numy / 896, 664, 448, 332, 224, 166 / 
data celle 6/122525-0 = 500 

data xydist / 3850.0 , 5350.0 , 3950.0 , 3950 / 


SLAT = 70. 

RE = 6378.273 
E2 = .006693883 
PI = 3.141592654 


E = sqrt(E2) 
C 
c Query for the SSM/I grid cell size. 
c 
print *, Enter the grid cell dimension:' 
oyetate es NAS) laser 
elstale yes 25.0) Ory 
Pritt. 5.50.0 kn 
read *,gtype 
C 
c Query for polar region of interest. 
c 


print *,’Enter the hemisphere of interest:' 
Pit oN Ont 

pints, 2 Sour, 

read *,ihem 
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c 
c Define the sign and meridian offset (delta) for the SSM/I grids. 
c 

if Ghem.eq.1) then 


SGN =1.0 
delta = 45. 
else 
SGN =-1.0 
delta = 0.0 
endif 
C 
c Query for translation type. 
C 
print *,’Enter one of the following transform functions:' 
print *, 1. Convert I,J to Latitude, Longitude’ 
print *,' 2. Convert Latitude, Longitude to I,J’ 
Teac. 1trans 
E 
c Start translation 
c 
if (itrans.eq.1) then 
€ 
c Obtain the I,J position of the grid cell to transform to Latitude 
c and Longitude 
c 
amt enterd.): 
read *,1,] 
c 


c Convert I,J pairs to x and y distances from origin. The grid will be 
c ‘flipped’ in the 'Y' direction. 


C 
x=((i-1)*cell(gtype))-(xydist(1,ihem)-cell(gtype)/2.) 
kk=numy(ihem,gtype)-G-1) 
y=((kk-1)*cell(gtype))-(xydist(2,ihem)-cell(gtype)/2.) 

C 


c Transform x and y distances to latitude and longitude 
c 
call mapxy (x,y,alat,alon,SLAT,SGN,E,RE) 


C 

c Transform radians to degrees. 

c 
alon=alon*180./PI 
alat=alat*180./PI 
alon=alon-delta 

C 

c Convert longitude to positive degrees 

c 


if (alon.le.0.0) alon=alon+360. 
if (alon.ge.360.0) alon=alon-360. 


SN 


c Print the latitude and longitude for the center of the I,J cell. 
¢ 
print *,alat,alon 
pause 
else 
e 
c Obtain the latitude and longitude pair and transform to cell where 
c that pair is located. 
E 
print *,’Enter latitude and longitude (positive values):' 
read *,lat,lon 


C 
c Transform degrees to radians 


i alat=abs(lat)*PI/180. 
alon=(lon+delta)*PI/180. 
: Transform latitude and longitude to x and y distances from origin 
° call mapll (x,y,alat,alon,SLAT,SGN,E,RE) 
: Convert x and y distances from origin to I,J pair (i,j) 
c 


li=nint((x+xydist(1 ,ihem)-cell(gtype)/2.)/cell(gtype))+1 
jjJ=nint((y+xydist(2,ihem)-cell(gtype)/2.)/cell(gtype))+1 


€ 
c Flip grid onentation in the 'Y’ direction 


Cc 

kk=numy(ihem,gtype)-(j- 1) 
Cc 
c Print the I,J location of the cell. 
c 

Oya bane ll de 

pause 

endif 
end 


SUBROUTINE MAPLL (X,Y,ALAT,ALONG,SLAT,SGN,E,RE) 


C$ DESCRIPTION: 


C$ This subroutine converts from geodetic latitude and longitude to Polar 

C$ Stereographic (X,Y) coordinates for the polar regions. The equations 

C$ are from Snyder, J. P., 1982, Map Projections Used by the U.S. 

C$ Geological Survey, Geological Survey Bulletin 1532, U.S. Government 
C$ Pnnting Office. See JPL Technical Memorandum 3349-85-101 for further 
C$ details. 
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C$. ARGUMENTS: 


CS 
C$ Variable Type I/O Description 
CS 


C$ ALAT #£REAL*4 I Geodetic Latitude (degrees, +90 to -90) 
C$ ALONG REAL*4 I Geodetic Longitude (degrees, 0 to 360) 
C$ xX REAL*4 O Polar Stereographic X Coordinate (km) 
Coe ¢ REAL*4 O Polar Stereographic Y Coordinate (km) 


C$ 

C$ Written by C. S. Morris - April 29, 1985 

C$ Revised by C. S. Morris - December 11, 1985 
C$ 

C$ Revised by V. J. Troisi - January 1990 

cS SGN - provides hemisphere dependency (+/- 1) 


REAL*4 X,Y,ALAT,ALONG,E,E?,CDR,PILSLAT,MC 


C$ 
5 DEFINITION OF CONSTANTS: 


C$ Conversion constant from degrees to radians = 57.29577951. 
CDR=57.2957795 1 
| AF shal ob 

C$  Pi=3.141592654. 
PI=3.141592654 

C Compute X and Y in grid coordinates. 
MetABS(ALAT).LT.PI/2.)GOTO 250 
X=0.0 
Y=0.0 
GOTO 999 

250 CONTINUE 
T=TAN(PI/4.-ALAT/2.)/((1.-E*SIN(ALAT))/(1.4+E*SIN(ALAT)))**(E/2.) 
IF (ABS(90.-SLAT).LT.1.E-5) THEN 
im O=2. “RE T/C +E) (14) 1-0-8) eZ.) 
ERSE 
Bees A PtP 130. 
eT AN(Pl/4 Si /(C1-B*SIN(SL))A 1 E*SINSE)))** (E72. ) 
mMe=COSSE)V SORT LO-E2*(GINGE) = 2)) 
MaO=RE*MC*T/TC 
END IF 
Y=-RHO*SGN*COS(SGN*ALONG) 
X= RHO*SGN*SIN(SGN*ALONG) 
999 CONTINUE 

END 


oo 


SUBROUTINE MAPXY (X,Y,ALAT,ALONG,SLAT,SGN,E,RE) 


(OG 4 ae a ae he ee ee ee ee He He Ae 2 he ie ie He eft He ee ie 2h Ae 2B FG ie RG Re oh 2 2 2 2 2 hc 2 2 a Ae oe oi oS ok Ae A AG os oe a te et 


CS DESCRIPTION: 


C$ This subroutine converts from Polar Stereographic (X,Y) coordinates 

C$ to geodetic latitude and longitude for the polar regions. The equations 

C$ are from Snyder, J. P., 1982, Map Projections Used by the U.S. 

C$ Geological Survey, Geological Survey Bulletin 1532, U.S. Government 
C$ Printing Office. See JPL Technical Memorandum 3349-85-101 for further 
C$ details. 


C$ ARGUMENTS: 

C$ Variable Type I/O Description 

CS xX REAL*4 I Polar Stereographic X Coordinate (km) 
Coan oe REAL*4 I Polar Stereographic Y Coordinate (km) 


C$ ALAT REAL *4 OQ Geodetic Latitude (degrees, +90 to -90) 
C$ ALONG REAL*4 OQ Geodetic Longitude (degrees, 0 to 360) 


CS 

ORS Written by C. S. Morris - April 29, 1985 

CS Revised by C. S. Morris - December 11, 1985 
C$ 

C$ Revised by V. J. Troisi - January 1990 

C$ SGN - provide hemisphere dependency (+/- 1) 
CS 


REAL*4 X,Y,ALAT,ALONG,E,E2,CDR,PI 


C$ 
C$ DEFINITION OF CONSTANTS: 


C$ Conversion constant from degrees to radians = 57.29577951. 
CDR=57- 205 7/251 
E2-b 
C$ Pi=3.141592654. 
PI=3.141592654 
I 
SL=SEAR AE! sO: 
200 RHO=SQRT(X**2+Y **2) 
IF (RHO.GT.0.1) GOTO 250 
ALAT=90.*SGN 
ALONG=0.0 
GOTO 999 
250 CM=COS(SL)/SQRT(1.0-E2*(SIN(SL)**2)) 
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T=TAN((PI/4.0)-(SL/(2.0)))/((1.0-E*SIN(SL))/ 
C(1.0+E*SIN(SL)))**(E/2.0) 
IF (ABS(SLAT-90.).LT.1.E-5) THEN 
T=RHO*SORT((1.+E)**(1.4E)*(1.-E)**(1.-E))/2./RE 
ELSE 
T=RHO*T/(RE*CM) 
END IF 
CHI=(PI/2.0)-2.0*ATAN(T) 
ALAT=CHI+((E2/2.0)+(5.0*E2**2.0/24.0)+(E2**3,0/12.0))*SIN(2*CHI)+ 
C((7.0*E2**2.0/48.0)+(29.0*E2**3/240.0))*SIN(4.0*CHI)+ 
C(7.0*E2**3.0/120.0)*SIN(6.0*CHI) 
ALAT=SGN*ALAT 
ALONG=ATAN2(SGN*#X,-SGN*Y) 
ALONG=SGN*ALONG 
999 CONTINUE 
END 
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program u, v calc 


2K 2k ok 2k ok ok ok 2k ok ok 2k ok ok 2k ok 2K 2k 2K ok OK 2K 2k OK ok ok 2K 2k Koi 2K 2K 2K 2 OK 2 2k 2K 2k 2K 2k 2K 2K ook ok 2c 2k ok ok ok 2K 2 ok oi ok aK ok ok ok ok ok 2K 2K 2k 2k 2K ok ok ok ok 


This program calculates the u and v components of motion from a list of 
user specified geodetic latitude and longitude postions. Results are output 
as an ASCII data file in units of cm/s. Appropriate lines in code must be 
modified to specify the filepath of input and output files. 


Variables: fname - input file name 
infile - input file path 
ofile - output file path 
lat,lon - latitude and longitude expressed as real numbers 
mlat - mean latitude between two positions 
u, V -uand v components of motion 
day - day of position (integer) 
length - length in days of series 
tmpd_- time between postions (in days) 


Mees ee ee 9G i eG ee I 


3K 2k ok 2k 2k ok ok ok ok 2k 2K 2k 2k ok 2k ok 2k 2k ok 2k 2k ok 2k 2k OK 2k ok 2K 2k ok 2 2 2K 2K 2k ok ok 2k 2 2K 2k OK 2K 2k ok ok 2k ok 2K 2 ok 2k 2 ok OK 2k ok OK ok ok 2 2k ok 2k ok ok ok ok ok 2 ok 2 


character* 10 fname 

character*70 ofile, infile 

real lat(60), lon(60), u(60), v(60),mlat 
integer day(60), length, tmpd 


c Construct input file path and open 


5 print *, ‘Enter input file name 
read (*,'(a)') fname 


print *, ‘Enter number of days in series’ 
read*, length 


print*, ‘Enter separation time period (days)’ 
read*, tmpd 


infile='Archs HardDisk:THESIS:Position Data:'//fname 
open (unit=15, file=infile,status='old',iostat=ier) 
if Ger.ne.Q) then 
print*,infile,'not accessible! <ret>' 
pause 
stop 
endif 


c Construct output file path and open 
ofile='Archs HardDisk:THESIS:Position Data:'//fname//’.uv' 
open (unit=18, file=ofile, status='new’, form='formatted’, 


+ 10Stat=ier) 
if (ier.ne.O)then 
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pant’, Pmor opening Zofile, <ret>' 
pause 
stop 
endif 
c Read in data 
read (15,10) (day), lat(1), lon), i=1, length) 
c Compute u, v components (km/time period) 
do 100 i=1, length-1 
yD) = (atl) = lat) 11 olla 
mlat = (lat(i) + lat(i+1))/2. 
call gcircl(mlat,lon@),mlat,lonG+1),uQ@)) 
if (latGi+1) .Jt. latG)) uG) = -uG) 


print®*, u(i), v@) 


100 continue 
pause 


c Convert u, v components to cm/s 
do 200 1=1, length-1 


m= ud) Clytimpa) (000 )*(100.)*(1./24,)5(1,/56000) 
v(i) = vG) * (1./tmpd)*(1000.)*(100.)*(1./24.)*(1./3600.) 


200 continue 
c Write to file 

write (18,20) ‘u’, 'v' 

write (18,20) ‘cm/s, cm/s 

write (18,30) (u@), vG), 1=1,length-1) 
KO etormat (10,3x,[5.2,2x,16.2) 
Poratonnat (1X, a; .5%.0) 


sy viommat (1x, 16.2... ,3x,f6.2) 


end 
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2 2k 2k 2 2k 2k 2k 2k 2k 2K oe KK ok ok ok ok ok ok OK ok ok kok ok ok ko ok kk ok ok ok ok ok ok Rk ok aK Kk KK a 
subroutine gcircl(lat] lon] ,lat2,lon2,gcircl) 
c This program computes great circle distance. 
c Inputs are latitude/longitude pairs in decimal degrees. 
c Output is circumferential distance between the points in kilometers. 
implicit real*8 (a-z) 
real*8 dpr 
real*4 latl, lonl, lat2, lon2, gcircl 
integer again 


dpr = 57.29577951308232 


c Compute difference in longitudes 
t = (lon2 - lonl) / dpr 
c Compute r 


dsinr = dsin(t) * dcos(lat2/dpr) 
r = dasin(dsinr) 


c Compute k 


dsink = dsin(lat2 / dpr) / dcos(r) 
k = dasin(dsink) 


c Compute angle D 


dcosd = dcos(r) * dcos(k - (latl / dpr)) 
d = dacos(dcosd) * dpr 


c Compute distance based on 111.6114 km/deg 
gcircl =d * 111.6114 
retum 


end 
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subroutine usetie (c , lhem ) 
BECECCeECECCCCCECECCCCCECCCEGECEEC CEC CEC CCE CECE ECECK CECE CCC CCCCECECCCECE 
C 
c USETIE - This subroutine computes coefficients for the ice algorithm 
C using the NASA-approved tiepoints. 
Inputs 


c  : Real array of coefficients. Refer to NODS User 
handbook for description of coefficients A-L. 


ihem : Integer denoting hemisphere (1=north, 2=south). Used 
in determining which set of coefficients to use. 
Internals 


TB...W: Open water tiepoints. One for each of the three channels 
used 1n the ice concentration algorithm (19V, 19H, 37V) 


TB...F: Firstyear tiepoints. One for each of the three channels 
used in the ice concentration algorithm (19V, 19H, 37V) 


TB...M: Multiyear tiepoints. One for each of the three channels 
used in the ice concentration algorithm (19V, 19H, 37V) 


DW _ : Difference between component tiepoints for open water. 
SW :Sum of the component tiepoints for open water. 

DF _ : Difference between component tiepoints for first year ice. 
SF : Sum of the component tiepoints for first year ice. 

DM_ : Difference between component tiepoints for multiyear ice. 


SM  : Sum of the component tiepoints for multiyear ice. 


Outputs 


c : Array of coefficients which have been computed from 
tiepoints. 


2 OO). 0 - OSS OQ OO (OO Oe GG OO) 4) OSG) _ Oe GOO Ga O21 eG 


ecececcecececcececececececececcceceecececececcececeeceeceecececececceCCcCCCCCS 


Fealc. 12) 
integer ihem 
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c Tiepoints for open water, first-year ice and multi-year ice 

c The first element of each tiepoint array represents the tiepoints 
c for the North Polar region; the second element is applied to the 
c South Polar region. 


c declare the two-element open water tiepoints 
real TBI9VW(2), TBISHW(2), TB37VW(2) 

c declare the two-element first-year tiepoints 
real TB19VF(2), TB19HF(2), TB37VF(2) 

c declare the two-element multi-year tiepoints 
real TB19VM(2), TB19HM(2), TB37VM(2) 


real DW(2) , DF(2) , DM(2) 
real SW(2) , SF(2) , SM) 


c First component of each tiepoint pair is for Northern Hemisphere 

c Second component of each tiepoint pair is of Southern Hemisphere 
c These tiepoints approved for use by the NASA Sea Ice Algorithm 
c Working Group (1/90) 


data TB19OVW / 175.3, 175.3/ 
data TBIOHW / 97.7, 97.7/ 
datasmelo V7 254.05 2512/ 
data TB19HF / 236.0, 241.7/ 
data TB19VM / 223.2, 223.2/ 
data TB19HM / 203.9, 203.9/ 
data TB37VW / 199.6, 199.6/ 
data TB37 VF / 250.0, 248.3/ 
data TB37VM / 186.3, 186.3/ 


DW(1) = TB19VW(ihem) - TB19HWGhem) 
DF(1) = TB19VFGhem) - TB19HF(ihem) 
DM(1) = TB19VM(ihem) - TB19HM(ihem) 
SW(1) = TB19VW(ihem) + TB19HW(ihem) 
SF(1) = TB19VF(ihem) + TB19HF(ihem) 
SM(1) = TB19VMGhem) + TB19HM(ihem) 
DW (2) = TB37 VW(ihem) - TB19V W(ihem) 
DF(2) = TB37VFGhem) - TB19VF(ihem) 
DM(2) = TB37VM(ihem) - TB19VM(ihem) 
SW(2) = TB37VW(ihem) + TB19VW(ihem) 
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SF(2) = TB37VFGhem) + TB19VF(ihem) 
SM(2) = TB37VM(ihem) + TB19VM(ihem) 


c Derive coefficients from tiepoints 


c(1) = DM(1)*DW(2) 
c(2) = DM(2)*SW(1) 


- DM(2)*DW(1) 
- DW(2)*SM(1) 
c(3) = DW(1)*SM(2) - 
c(4) = SM(1)*SW(2) - 
c(5) = DF(1)*(DM(2) - 
1+ = DM(1)*(DW(2) - 
c(6) = DF(2)*(SM(1) - 


DM(1)*SW(2) 

SM(2)*SW(1) 

DW(2)) + DW(1)*(DF(2) - DM(2)) 
DF(2)) 

SW(1)) + DW(2)*(SF(1) - SM(1)) 


eZ) CSW4 1) = SG) 


c(7) = DF(1)*(SW(2) - 


SM(2)) + DW(1)*(SM(2) - SF(2)) 


fs DMC) *(SFR)-Sw)) 


c(8) = SF(2)*(SW(1) - 


SM(1)) + SW(2)*(SM(1) - SF(1)) 


1+ SM(2)*(SF(1) - SW(1)) 

c(9) = DF(2)*DW(1) - DF(1)*DW(2) 
c(10) = DW(2)*SF(1) - DF(2)*SW(1) 
c(11) = SW(2)*DF(1) - DW(1)*SF(2) 
c(12) = SF(2)*SW(1) - SF(1)*SW(2) 


return 
end 
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subroutine 1ce25 ( tb,c,maxlen,idim1,idim2,ice,mask,icepar ) 
CCCCCECCCCCCCCCCCECCCECCECECECEECCCCCCCCCCCOBECECCEGEBEECCCCEEe cCCCCCU CE 
c ICE25 - Compute ice concentrations from NSIDC CD-ROM data at 25 km 
resolution. 


Q 


Software 1s modified from the Vince Troisi's original 50 km version by 
W.A. Wright, Naval Postgraduate School. 


Algonthm 
The ice concentration algorithm is an adaptation of the 
algorithm designed by the NODS Software Group (refer to the 
NODS User Handbook, Appendix CC, pp 22-26). The algorithm 
uses the latest SSM/I ttepoints which were determined by Don 
Cavalier and Koni Steffen. 


A weather filter (GR > .05) is applied to correct for false 
concentrations computed in regions of open water. 


Computed concentrations greater than 120% and less than -20% 
are considered suspect and the corresponding cell 1s flagged 
as missing (-99). 

Inputs 


tb : buffer of TB cell values. This buffer includes all 
SSM/I TB parameters excluding the 85GHz channels. 


idim1 : The number of parameters in a SSM/I cell (ie., 5) 

idim2 : The number of cells in a SSM/I Tb grid line. 

maxlen: The maximum length of tb 

mask : A line from the landmask for 25 km grid. The mask 
line corresponds to the ice grid line to be 
computed. (1=land; O=water) 

icepar: The ice concentration type to compute (1=total, 
2=multiyear) 

Internals 


c : Array of coefficients used to compute ice concentrations. 
Refer to NODS Handbook Volue II. Appendix CC pp 23-24. 


ioff : Offset pointers used to point to Tb values used in 
this ice algorithm (19GHz V, 19GHz H, 37GHz V). 


ee Ga ee ey on) CoG) OC) GC). 1) 4G) On OO) ) Civ GQ Qn. - OO ONG GeO. © ©) Gi QQ: G89 aA O80 O°. aa 6 


ipoOs : position in tb array which 1s to be summed 
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sum : Real register used for summing four adjacent Tb values 


count : The number of Tb values not flagged as missing. A Tb value 
of zero (0) is a missing Tb value. 


fyr : Integer register which stores the calculated first-year 
ice concentration 


myr_: Integer register which stores the calculated multi-year 
ice concentration 


tot : Integer register which stores the calculated total ice 
concentration 


fyrb : Two - byte buffer equivalenced to fyr register. Only the 
first byte is the significant part of the register. 


myrb : Two - byte buffer equivalenced to myr register. Only the 
first byte is the significant part of the register. 


GR __ : Gradient ratio 
PR _ : Polarization ratio 
Outputs 
ice : Line of 25 km ice concentration grid. Values are 
represented as a single byte percentage ranging 
from 0 - 100%. Cells over land flagged with -88; 


cells for which concentrations were suspect or 
not computed are flagged with -99. 


COC ae) GO OG) ON Od He Ce oC pCa eC a Cr Ga a Os a re CN Gt Ge ©) eG 


Gece CEEECECECCCECCCCECCCECC CEC ECC CCCCECECC ECC CCC ECC CEE CCECCCECECCCCe CCC 
integer idim1,idim2 
integer maxlen 
integer*2 tb(maxlen*idim1) 


c Offset from start of cell to three parameters used in ice 
c calculations for two lines of data 


imeger*2 ioff(3) 

integer*2 fyr,myr,tot 

integer*1 fyrb(2),myrb(2) 

integer*1 ice(idim2),mask(idim2),fyrice,myrice 


logical GOOD 


c Three TB parameters for ice concentrations 
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real BT19V,BT19H,BT37V 
real CCl) 

real GR,PR,PRGR 

real ANF,ANM,DD 


equivalence (fyr, fyrb(1)) 
equivalence (myr, myrb(1)) 


c Offset pointers for Brightness temperature cells used in algorithm 
data ioff /0,1,3/ 
c Begin processing buffer 


ii =O 
do 1000 1 = 1,idim1 *idim2,idim1 


c Increment pointer for next Ice cell. 
i=) 
c Calculate Ice Concentration _ 
GOOD = .true. 
c Set ice pixel to land(-88) if over land 
if (mask(ii).eq.1) then 
ice(11)=-88 
else 
c Determine Ice Concentration of the ice pixel 
do 80 l=1,3 
ipos=i+10ff(1) 
tbval = (float(tb(ipos)))/10. 
if(l.eq.1) BT19V = tbval 
ied. 2) Bi 19H = tbval 


if(l.eq.3) BT37V = tbval 
SO continue 
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c Set ice cell to missing if any TB parameters are missing 


BEY BOO 0 or B19 EO 0.0 <or 


si BT37V.EQ.0.0 ) then 
fyrice =-99 
myrice =-99 
ice(iil) =-99 
else 


c Compute open water weather filter 
GR = (BT37V - BT19V)/((BT37V + BT19V) 


c Set ice cell to 0% concentration if weather filter > .05 
c Filter value of 0.05 may be altered at a future date! 


if (GReet..05) then 
fycice. — 0 
msmice — 0 
ice = 0 
else 
c Apply noniterative ice algorithm 
c Compute polarization ration PR 
PR =(BTI19V - BT19H)/(BT19V + BT19H) 
PRGR = PR*GR 
ANF = c(1) + c(2)*PR + c(3)*GR + c(4)*PRGR 
ANM =c(9) +c(10)*PR + c(11)*GR + c(12)*PRGR 
DD =c(5) + c(6)*PR + c(7)*GR + c(8)*PRGR 
c Compute First-Year, Multi-Year, and Total Ice Concentrations 


fyr = 1f1x(100.0*(ANF/DD)) 
myr = 1f1x(100.0*(ANM/DD)) 
tot = fyr + myr 


c Filter out values outside the range of -20 to 120 
c This filter may be altered at a future date 


if (fyr .1t. -20 .or. 


et myr .|t. -20 .or. 

tot .lt. -20 .or. 

- fyr .gt. 120 .or. 

= myr .gt. 120 .or. 

: tot .gt. 120) then 
fyrice = -99 
myrice =-99 
ice(il) = -99 
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GOOD = fase. 
endif 


c Place concentration values into byte data representations 


if (GOOD) then 
fyrice = fyrb(2) 
myrice = myrb(2) 
if GCepatequriceal) = fyrice + mymice 
if (acepar.eq.2) ice(ii) = myrice 


c Set concentration to 0 if computed concentration between -20 and 0 
c Set concentration to 100 1f computed concentration between 100 and 120 


ieace(mmlt 0) 1ceGi) =O 
if (ice(ii) .gt. 100) iceGi) = 100 

endif 

GOOD =.true. 

endif 
endif 
endif 
1000 continue 


c Return computed ice concentration array, ice 


return 
end 


Viz 
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